Causal inference method and system for agricultural meteorological disaster loss assessment

Abstract

The invention discloses a causal inference method and system for agricultural meteorological disaster loss assessment, and belongs to the field of agricultural meteorological disaster prevention and reduction. The method comprises the following steps: acquiring multi-source spatio-temporal data including time sequence data, static data, disaster data and crop yield observation values; performing time sequence processing on the time sequence data to generate a depth time sequence feature, and generating an attention weight corresponding to each time step in the time sequence data through a multi-head self-attention mechanism; the depth time sequence features and the static data are input into a first regression device for training, time step weighting is performed on the time sequence data according to the attention weight, the time sequence data and the static data are input into a second regression device for training, an SHAP interpreter is adopted to interpret and analyze the second regression device, and the influence of the environment features on the yield is quantified; the disaster data is used as an intervention variable, the time sequence data and the static data are used as covariables, and the causal effect evaluation is performed according to the intervention variable, the covariables and the crop yield observation value to obtain corresponding crop yield variations under different disaster levels.

Description

Technical Field

[0001] This invention relates to the field of agricultural meteorological disaster prevention and mitigation, specifically to a causal inference method and system for assessing agricultural meteorological disaster losses. Background Technology

[0002] Against the backdrop of intensifying global climate change, agricultural meteorological disasters such as droughts, floods, and frosts are occurring frequently, posing a severe challenge to food security. Timely and accurate assessment of the specific losses to crop yields caused by these disasters is crucial for developing disaster mitigation measures, processing agricultural insurance claims, and optimizing disaster response policies. To quantify the losses caused by a specific disaster in a given year (e.g., 2023), it is necessary to compare the observed yield under actual disaster conditions with the hypothetical yield under conditions where no disaster occurred. However, in actual production, once a disaster occurs, the counterfactual state of no disaster cannot be directly observed, which constitutes a fundamental technical challenge in disaster loss assessment.

[0003] Currently, traditional crop yield prediction methods are mainly divided into the following two categories, but both have significant limitations: 1. Relying on crop mechanistic models such as WOFOST, these methods predict potential yields by simulating the interaction between environmental factors such as light, temperature, water, and soil and crop physiological processes. Although crop mechanistic models have transparent and clear mechanisms and strong interpretability, they are not suitable for regional scales due to the difficulty in obtaining model parameters and data.

[0004] 2. Relying on purely data-driven models, this type of method uses historical weather, remote sensing, and crop yield data to build predictive models, avoiding complex mechanistic modeling. However, these deep learning methods are usually "black box" structures, unable to trace the intermediate process from input values ​​to predicted values, and cannot easily explain the relationship between environment and yield.

[0005] Therefore, it is crucial to accurately assess the causal impact of agricultural meteorological disasters on crop yields and to develop targeted disaster mitigation measures. Summary of the Invention

[0006] In view of this, the main objective of the present invention is to provide a causal inference method and system for assessing agricultural meteorological disaster losses, in order to at least partially solve the above-mentioned technical problems.

[0007] To achieve the above objectives, as the first aspect of this application, a causal inference method for assessing agricultural meteorological disaster losses is proposed, comprising the following steps: S1: Acquire multi-source spatiotemporal data, which includes time-series data, static data, disaster data, and crop yield observations. The time-series data includes meteorological and remote sensing data. S2: Perform time-series processing on the time-series data to generate deep time-series features, and generate attention weights corresponding to each time step in the time-series data through a multi-head self-attention mechanism; S3: Input the deep temporal features and the static data into the first regressor for training. At the same time, weight the temporal data according to the attention weight and input it into the second regressor along with the static data for training. Use the SHAP interpreter to interpret and analyze the second regressor to quantify the impact of environmental features on yield. S4: The disaster data is used as the intervention variable, and the time series data and the static data are used as covariates. The causal effect is evaluated based on the intervention variable, covariates and crop yield observations to obtain the crop yield change corresponding to different disaster levels.

[0008] In one possible implementation, the time-step weighting of the time-series data based on attention weights described in S3 specifically involves: multiplying the attention weight corresponding to each time step with the observation values ​​of all time-series data within that time step to generate an attention-weighted time-series feature matrix; standardizing the static data to obtain static features, expanding them along the time dimension, and concatenating them with the attention-weighted time-series feature matrix to obtain the attention-weighted raw input; and flattening the attention-weighted raw input along the time-step dimension to form a two-dimensional feature matrix for training the second regressor.

[0009] In one possible implementation, a leave-one-year cross-validation strategy is used to train the first and second regressors, and the average of the evaluation metrics from the n training iterations is taken as the final evaluation of the model performance; where n is the total number of years. By selecting data from one year as the validation set and data from the (n-1)th year as the training set, n independent training iterations are performed to train the first and second regressors, and the average of the evaluation metrics from the n training iterations is taken as the final evaluation of the model performance.

[0010] In one possible approach, a causal forest model is used to assess causal effects and obtain individual treatment effect values ​​under different disaster levels, whereby the individual treatment effect values ​​represent the changes in crop yield corresponding to different disaster levels.

[0011] In one possible implementation, the causal forest model is trained using a leave-one-year cross-validation strategy before step S4, during which the validation set corresponding to disaster types that do not appear in the training set is filtered out.

[0012] In one possible implementation, S4 is followed by a visualization step: plotting the average causal effect curve and the spatial distribution map of individual treatment effects based on the causal effect assessment results.

[0013] In one possible implementation, the causal forest model is trained using a leave-one-year cross-validation strategy before step S4, during which the validation set corresponding to disaster types that do not appear in the training set is filtered out.

[0014] In one possible implementation, the multi-source spatiotemporal data is preprocessed before S2. The preprocessing includes the following steps: reshaping the meteorological and remote sensing data to form a two-dimensional matrix, the two-dimensional matrix being set as the number of samples × the time step, and standardizing the two-dimensional matrix to obtain preprocessed time-series data.

[0015] In one possible implementation, both the first regressor and the second regressor are XGBoost regression models.

[0016] As a second aspect of this application, a causal inference system for assessing agricultural meteorological disaster losses is also proposed, the system comprising: The data acquisition unit is used to acquire multi-source spatiotemporal data, which includes time-series data, static data, disaster data, and crop yield observations. The time-series data includes meteorological and remote sensing data. The feature extraction unit, connected to the data acquisition unit, is used to perform time-series processing on the time-series data to generate deep time-series features, and to generate attention weights corresponding to each time step in the time-series data through a multi-head self-attention mechanism. The dual-path regression unit includes a first regressor and a second regressor, which are connected to the feature extraction unit. The unit is used to input the deep temporal features and the static data into the first regressor for training, and at the same time, to weight the temporal data according to the attention weights and input it into the second regressor along with the static data for training. The interpretation unit uses the SHAP interpreter to interpret and analyze the second regressor, quantifying the impact of environmental characteristics on yield. The causal assessment unit uses the disaster data as an intervention variable. It uses the time-series data and the static data as covariates, and combines the intervention variable and crop yield observations to assess the causal effect and output individual treatment effect values. The individual treatment effect values ​​represent the changes in crop yield corresponding to different disaster levels.

[0017] One possible implementation includes: a preprocessing unit for reshaping the meteorological and remote sensing data to form a two-dimensional matrix, the two-dimensional matrix being set as the number of samples × time step, and for standardizing the two-dimensional matrix to obtain preprocessed time-series data.

[0018] As a third aspect of this application, an electronic device is also proposed, including a processor and a memory, the memory storing a computer program, characterized in that the processor executes the computer program to implement the steps of the method described above.

[0019] Based on the above technical solution, it can be seen that the causal inference method and system for assessing agricultural meteorological disaster losses proposed in this application has at least one of the following beneficial effects compared with the prior art: 1. By constructing a dual-path regression prediction system, on the one hand, time-series data is processed to generate deep time-series features, ensuring the model's ability to learn complex nonlinear time-series relationships and its prediction accuracy; on the other hand, attention weights are used to weight the original time-series data at the time step level, generating explicit and interpretable features for prediction. This dual-path design not only improves the robustness of prediction through comparative verification, but its second path also directly traces the model's decision-making mechanism back to specific original environmental variables, solving the black-box problem of deep learning models in agricultural applications, and making the final yield prediction highly reliable and interpretable.

[0020] 2. By performing extensive training on the time-series deep learning model to stabilize features, and then using leave-one-year cross-validation to evaluate model generalization, robustness and repeatability under inter-year data distribution shifts are ensured. The leave-one-year cross-validation method is used to evaluate the model's cross-year generalization performance, effectively solving the problem of data distribution changes caused by inter-year climate fluctuations and adjustments to planting systems. This ensures that the model can maintain stable and reliable performance in disaster assessment in unknown years, improving the consistency and reliability of the technology in agricultural applications on a large spatiotemporal scale.

[0021] 3. Employing a causal forest model, using time-series and static data as covariates and disaster data as the intervention variable, combined with crop yield observations for causal effect assessment, effectively isolates confounding environmental influences, revealing the direct crop yield changes caused by different disaster levels. Simultaneously, the SHAP interpreter is used to interpret and analyze the second regressor, quantifying various environmental characteristics. Ultimately, the causal forest model provides reliable causal inferences for disaster loss quantification, while SHAP analysis offers interpretable insights into the impact of environmental factors. These two approaches complement and support each other, together forming a complete and interpretable agricultural meteorological disaster loss assessment and analysis system, from causal quantification of disaster damage to detailed analysis of environmental contributions. Attached Figure Description

[0022] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is a flowchart illustrating a causal inference method for assessing agricultural meteorological disaster losses according to the present invention; Figure 2 This is a schematic diagram of an electronic device according to the present invention. Detailed Implementation

[0024] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to specific embodiments and accompanying drawings.

[0025] The terminology used in this invention is for the purpose of describing particular embodiments only and is not intended to limit the embodiments of the invention. The singular forms “a,” “the,” “and,” and “the” as used in the embodiments of the invention and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise.

[0026] Existing methods for assessing agricultural meteorological disaster losses suffer from two main shortcomings: first, they cannot accurately pinpoint the key impacts identified by the model to specific crop growth stages and their corresponding original environmental variables; second, in actual production, the counterfactual state of no disaster cannot be directly observed once a disaster occurs, leading to an inability to obtain the true yield loss. This constitutes a fundamental technical challenge in disaster loss assessment. To address these issues, in-depth research has revealed that training two regressors in parallel, synergistically utilizing deep time-series features and interpretable weighted time-series features, ensures the accuracy and robustness of yield predictions.

[0027] The present invention will be further illustrated below through specific embodiments. It should be noted that the following embodiments are merely illustrative and not intended to limit the scope of the invention. All other embodiments obtained by those skilled in the art based on the embodiments shown below without inventive effort are within the scope of protection of the embodiments of the present invention.

[0028] Taking the assessment of agricultural meteorological disaster losses in rice-growing areas of Northeast China as an example, this paper describes in detail the specific implementation process of the entire technical solution.

[0029] Therefore, as Figure 1 As shown, the inventors have proposed a causal inference method for assessing agricultural meteorological disaster losses, comprising the following steps: S1: Acquire multi-source spatiotemporal data. In this embodiment of the application, multi-source spatiotemporal data for 11 years from 2010 to 2020 in Northeast China (including Heilongjiang, Jilin and Liaoning provinces) are collected. These data come from meteorological stations, remote sensing satellites and agricultural statistics departments.

[0030] Multi-source spatiotemporal data includes time-series data, static data, disaster data, and crop yield observations. The time-series data is obtained from 8-day synthetic data provided by the MODIS remote sensing satellite. The original data structure of the time-series data is a three-dimensional tensor [number of samples × time step × number of features] = [200 × 29 × 8], covering 8 days of synthetic data for the annual rice growing season (days 65 to 297, a total of 29 time steps, 200 county-level samples). The time-series data includes meteorological data and remote sensing data; meteorological data includes daily average temperature, cumulative precipitation, solar radiation, relative humidity, and wind speed; remote sensing data includes Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Solar Induced Chlorophyll Fluorescence (SIF). Static data includes soil type data, soil organic matter content data, soil pH data, topographic slope data, and elevation data. These characteristics remain largely unchanged throughout the rice growing season. Crop yield observations are obtained by collecting rice yield data per unit area from each sample, in kilograms per hectare. In addition, the multi-source spatiotemporal data also includes geographic location data, which records the geographic coordinates of each sample, including latitude and longitude. For the latitude and longitude vectors, two one-dimensional arrays named "latitude" and "longitude" are extracted. The latitude vector stores the latitude coordinates of each sample, and the longitude vector stores the corresponding longitude coordinates, both represented in decimal degrees, for geographic location verification to ensure that the coordinates of all samples are located within the study area.

[0031] In this embodiment of the application, all the above data are unified onto the same geographic grid, and the required feature variables, all features and labels of all sample points in the same year are packaged and stored as an NPZ file. The NPZ file is a very common and standard data storage file format in the field of machine learning, used to save compressed archive files of multiple multidimensional arrays.

[0032] In one possible implementation, the multi-source spatiotemporal data is preprocessed before S2. The preprocessing includes the following steps: reshaping the meteorological and remote sensing data to form a two-dimensional matrix, the two-dimensional matrix being set as the number of samples × time step, and performing mean imputation and standardization on the two-dimensional matrix to obtain preprocessed time-series data; and standardizing the static data to obtain static features.

[0033] In this embodiment, since preprocessing algorithms such as missing value imputation and standardization typically require input of two-dimensional tabular data (rows = sample units, columns = features), the original three-dimensional time series data needs to be reshaped first. Therefore, the three-dimensional structure of sample number × time step × feature number is converted into a two-dimensional structure of (sample number × time step) × feature number, that is, the 29 time steps of each sample are split into independent sample-time step units. The original three-dimensional tensor [200×29×8] is merged along the sample number and time step dimensions to calculate a new row dimension of 200×29=5800 (a total of 5800 sample-time step units), while the column dimension remains unchanged with 8 features, finally generating a two-dimensional matrix [5800×8]. This integrates the time series data scattered at different time steps into a unified two-dimensional table, ensuring that subsequent preprocessing algorithms can accurately imput the missing values ​​of each sample-time step.

[0034] In this embodiment, since static data has no time dimension, there is no need for reshaping or dimension recovery. Only standardization is needed to eliminate dimensional differences. The specific operation process is as follows: The original static data is a two-dimensional matrix [200×5] (200 samples, 5 static features). Using the same standardization tool as for time-series data (StandardScaler), mean-variance standardization is performed separately for each static feature column, such as soil organic matter content in g / kg and altitude in m, ensuring that the mean of each static feature is 0 and the variance is 1. This is to avoid excessive focus on high-value features and neglect of low-value but important features (such as the effect of pH on rice root absorption) due to large differences in altitude (range 100-1000m) and soil pH (range 5.5-8.5). The resulting preprocessed static features are still a two-dimensional matrix [200×5], which will be used as subsequent data input. Furthermore, the data quantity of each sample is checked and ensured to be consistent, and the spatiotemporal correspondence is correct. If inconsistencies are found, an error log is recorded and the processing flow is terminated. If they match, all the parsed data components are packaged and returned to provide structured input data for subsequent data preprocessing steps. This process ensures the effective integration of multi-source heterogeneous data and establishes a reliable data foundation for subsequent feature extraction and causal inference.

[0035] S2: Perform time-series processing on the time-series data to generate deep time-series features, and generate attention weights corresponding to each time step in the time-series data through a multi-head self-attention mechanism; In this embodiment, preprocessed temporal data is input into a trained BiLSTM-Attention model. The model's network structure is configured as a bidirectional LSTM layer: two layers, each with 64 hidden units, capturing temporal dependencies from both forward and backward directions. The multi-head self-attention layer includes eight attention heads, each independently calculating attention weights. Residual connections: input and output are fused through a fully connected layer, using ReLU activation, and the hidden state of the last layer is extracted as a deep temporal feature with a shape of [200, 128]. The attention weight matrix [200, 29, 8] is output from the multi-head self-attention layer. Here, the eight features [200, 29, 8] can be understood as eight attention heads, each calculating a weight for 29 time steps (assuming aggregation along the feature dimension, or defaulting to each head outputting a weight for each time step). The mean is calculated along the head dimension, and the opinions of the eight heads are combined to obtain the final consensus weight for each time step. We get [200, 29, 1], which is the single attention weight for each sample and each time step, with the shape [200, 29] (the last dimension is usually compressed).

[0036] In one possible implementation, as described in S3, the temporal data is weighted according to attention weights at each time step by aggregating (e.g., averaging) the weights of the multi-head outputs in the feature dimension and / or head dimension to generate an attention weight vector corresponding to each time step, which has the shape of [number of samples, number of time steps] or [number of samples, number of time steps, 1]. Subsequently, this weight vector is used to perform time-step weighting on the original temporal feature matrix. Specifically, the attention weight corresponding to each time step is multiplied by the observation values ​​of all temporal data within that time step to generate an attention-weighted temporal feature matrix. The shape of the attention-weighted temporal feature matrix is ​​[B, T, F], where B represents the number of samples, T represents the time step, and F represents the number of temporal features. The shape of the static features of the multi-source spatiotemporal data is [B, S], where S represents the number of static features. The static features are expanded in the time dimension to obtain the expanded static features [B, T, S]. The expanded static features are concatenated with the attention-weighted temporal feature matrix to obtain the attention-weighted original input of shape [B, T, F+S]. The features of each of the T time steps of the multi-source spatiotemporal data are flattened in chronological order into a one-dimensional feature vector of length T×(F+S), thus forming a two-dimensional matrix of shape [B, T×(F+S)] for training the second regressor.

[0037] In one possible implementation, both the first and second regressors are XGBoost regression models.

[0038] In this embodiment, when the attention weights are multiplied element-wise with the temporal data [200, 29, 4], the actual result is that the weight vector [200, 29, 1] × the temporal data [200, 29, 4] = the attention-weighted temporal feature matrix [200, 29, 4]. One weight at each time step is broadcast and applied to the four features at that time step; this is the meaning of time-step weighting. The original static feature matrix [200, 3] is copied and expanded along the time dimension. Specifically, the three static feature values ​​of each sample are copied 29 times along the time dimension to form the expanded static feature matrix [200, 29, 3]. The expanded static feature matrix [200, 29, 3] and the weighted temporal feature matrix [200, 29, 4] are concatenated along the feature dimension to form the complete attention-weighted original input [200, 29, 7]. After concatenation, the feature vector at each time step consists of 4 weighted temporal variables and 3 static variables arranged sequentially, for a total of 4 + 3 = 7 features. Therefore, the final complete attention-weighted original input is a three-dimensional tensor with the shape [200, 29, 7]. After concatenation, the feature vector at each time step contains 7 elements: the first 4 are weighted temporal features, and the last 3 are extended static features. The three-dimensional attention-weighted original input [200, 29, 7] is flattened into a two-dimensional matrix. Specifically, for each sample (200 samples), its 29 time-step features are arranged sequentially to form a one-dimensional feature vector. The length of the flattened feature vector for each sample is T × (F + S) = 29 × 7 = 203. This transformation converts the spatiotemporally structured input into a traditional two-dimensional tabular data format (sample × feature), allowing the XGBoost regressor to process it directly. Using the flattened 2D matrix [200, 203] as the feature input and the known crop yield of the corresponding sample as the label, a second XGBoost regressor (i.e., the second-path regressor) is trained. By directly applying attention weights to the original features, the model's focus on key time steps and key factors can be traced back to the specific original variables, solving the black box problem of deep learning models. The flattening operation fully preserves the temporal sequence information, enabling the regressor to learn patterns in the time dimension, rather than just the statistical regularities of feature values.

[0039] In this embodiment, taking 2020 as the validation year as an example, the training set consists of 10 years of data from 2010 to 2019, with a sample size of "10 × 200 = 2000 samples". Two input features for the regressors need to be generated: ① Input for the first regressor training: For the 2000 samples of time-series data, the pre-trained BiLSTM-Attention model extracts deep time-series features [2000 × 128], which are then concatenated with static features [2000 × 3] to form [2000 × 131]; ② Input for the second regressor training: For the 2000 samples, an attention-weighted two-dimensional feature matrix [2000 × 203] is generated according to step S3, with attention weights generated synchronously by the BiLSTM-Attention model. The validation set consists of 2020 data, with a sample size of 200. Similarly, two input features for the regressors need to be generated, with the same format as the training set, only differing in sample size, for subsequent model prediction and performance evaluation. Next, the first regressor is loaded with the following parameters: n_estimators=100, learning_rate=0.1, max_depth=6. Using the training set input [2000×131] as features and the corresponding "actual rice yield observation value" (unit: kg / hm²) as labels, the model is trained using a function. Training stops when the validation set loss shows no decrease for 10 consecutive rounds to avoid overfitting. After training, the first yield prediction model is saved as model file A for subsequent yield prediction of samples. The second yield prediction model is then trained, with initialization parameters identical to the first regressor to ensure no parameter bias. Using the training set input [2000×203] as features and the same "actual rice yield observation value" as labels, the model is trained using a function. After training, the second yield prediction model is saved as model file B, thus completing model training. The SHAP interpreter is used to perform feature importance analysis on the second regressor to explain the contribution of each feature to yield prediction.

[0040] In this embodiment, an interpreter based on SHAP theory, specifically TreeExplainer for tree models, is used to perform importance analysis on the XGBoost regressor of the second regression path. The specific implementation process is as follows: 1. SHAP Interpreter Initialization and Configuration First, initialize the TreeExplainer interpreter and bind it to the trained XGBoost model. Specific steps include: configuring the interpreter to load the best-performing XGBoost model trained using the second regression path; creating a TreeExplainer instance, specifying model parameters and calculation methods; and setting the calculation method for feature expectation to conditional expectation based on the training data. For the background dataset selection, 100 representative samples are randomly selected from the training set as the explanatory background dataset. Selection criteria include: covering samples with different yield levels (high, medium, low yield), including samples with different cold disaster levels, ensuring representative geographical distribution (plains, hills, mountains), and maintaining the integrity of feature distribution.

[0041] 2. SHAP value calculation process Using all samples from the validation set (100 county-level samples) as the dataset to be explained, SHAP values ​​are calculated as follows: For each validation sample, its SHAP value is calculated across all feature dimensions; parallel computing is used to optimize processing speed, and the SHAP values ​​of all samples are calculated at once.

[0042] Based on the SHAP values ​​of all samples, they are arranged into a SHAP value matrix. The rows of this matrix represent 100 validation samples, and the columns represent 203 feature dimensions. Each element represents the contribution of a feature to the prediction result of a sample.

[0043] 3. Feature Name Generation and Mapping To facilitate interpretation, a clear naming system needs to be established for the 203-dimensional features, corresponding to their agronomic meanings: Based on the structure of the attention-weighted raw input, a systematic feature naming system is established: Time series feature naming: "variable_name_t_timestep", example: "NDVI_t18" represents the normalized vegetation index at time step 18; "Temperature_t19" represents the temperature at time step 19; "Precipitation_t15" indicates the precipitation at the 15th time step; "SIF_t20" indicates the solar-induced chlorophyll fluorescence at time step 20.

[0044] Static feature naming: "variable_name_static", example: "OrganicMatter_static" indicates the soil organic matter content; "pH_static" indicates the soil pH level; "Elevation_static" represents the altitude.

[0045] This results in a complete feature name mapping table, through which abstract multidimensional features can be transformed into variable names with clear agronomic meaning.

[0046] For a given feature, first take the absolute value of the feature's SHAP value across all samples, then calculate the average of these absolute values. This value represents the average magnitude of the feature's influence on the model's prediction results; it is a non-negative scalar. The features are ranked in descending order of importance score to obtain the feature importance ranking.

[0047] 4. Implementation Example Analysis Results: The top 10 features by importance are: SIF_t20: 0.148 (chlorophyll fluorescence at heading stage); Temperature_t19: 0.132 (temperature during heading stage); NDVI_t16: 0.115 (vegetation index at jointing stage); SIF_t18: 0.102 (chlorophyll fluorescence at the early heading stage); Temperature_t21: 0.095 (temperature during the later stages of heading); Precipitation_t15: 0.088 (precipitation during the jointing stage); OrganicMatter_static: 0.075 (soil organic matter); NDVI_t12: 0.068 (late tillering vegetation index); pH_static: 0.054 (soil pH); Elevation_static: 0.048 (elevation); 5. Co-validation with attention weights To verify whether the key reproductive period focused by the attention mechanism is consistent with the feature importance, the SHAP feature importance and attention weight are compared and analyzed. The average attention weight for each time step is calculated, and the sum of the SHAP values ​​for the corresponding features at each time step is calculated. In this embodiment, the 29 time steps are divided into several reproductive stages (such as tillering stage, jointing stage, and heading stage), and the average attention weight of all time steps within each stage is calculated.

[0048] Analysis showed that the model exhibited the highest attention level (mean attention weight 0.185) and the largest characteristic influence (mean SHAP contribution 0.126) during the heading stage (time steps 18-22); followed by the jointing stage (time steps 12-17) (mean attention weight 0.124, mean SHAP contribution 0.089); and relatively lower attention level during the tillering stage (time steps 8-11) (mean attention weight 0.087, mean SHAP contribution 0.062).

[0049] The analysis results are highly consistent with rice physiology. The heading stage is the critical period for yield formation and is most sensitive to environmental stress. Chlorophyll fluorescence (SIF) directly reflects photosynthetic efficiency. Temperature within a suitable range promotes growth, while exceeding the threshold causes stress.

[0050] Under different agricultural meteorological disaster scenarios, the attention weights in this invention can adaptively establish strong correlations with the key crop growth cycles most sensitive to the corresponding disaster stress. For example, in the assessment of low-temperature chilling injury in rice, the attention weights will significantly focus on the heading stage (such as time steps 18-22), which is sensitive to low temperatures, because this stage directly affects floret fertility and seed setting rate. In drought assessment, the model may assign higher weights to the crop's critical water period (the critical water period refers to the key growth stage during crop growth and development that is most sensitive to water deficiency, and once water is lacking, it will lead to a significant decrease in yield that is difficult to compensate for by later watering), such as the tasseling and silking stage of maize or the flowering and pod-setting stage of soybean. This dynamic correlation stems from the model learning the correspondence between specific disaster types and specific stress critical periods from historical data containing different interventions during the training process. The model not only provides a unified time series analysis framework, but also automatically and accurately identifies and traces back to the core reproductive stage that caused the damage based on specific disaster interventions. This generates causal attributions that are highly targeted and agronomically interpretable, providing precise time windows and key variable targets for the formulation of differentiated disaster reduction strategies for different disasters.

[0051] In one possible implementation, a leave-one-year cross-validation strategy is used to train the first and second regressors. Data from one year is selected as the validation set, and data from n-1 years are used as the training set. A total of n independent training iterations are performed to train the first and second regressors. The average of the evaluation metrics from these n training iterations is taken as the final evaluation of the model performance, where n is the total number of years. Ultimately, by comparing the comprehensive performance of each model on the validation sets of different years, the aim is to determine the optimal model configuration for each candidate training set and fairly select the predictive model with better cross-year generalization ability.

[0052] In this embodiment, the model is iterated 11 times annually (2010-2020). In each iteration, data from a single year is used as the validation set, and data from the remaining 10 years are used as the training set. The model performance is evaluated using two metrics: the coefficient of determination (R²) and the root mean square error (RMSE). The coefficient of determination (R²) represents the degree to which the independent variable explains the variation in the dependent variable, and the root mean square error (RMSE) is the ratio of the square root of the number of sample points to the sum of the squares of the deviations between the predicted and actual values. , Where N represents the sample size. This represents the actual value of sample i. R² represents the model's predicted value for sample i. The closer R² is to 1 and the smaller the RMSE, the better the model's performance. The performance of the first regressor (R²=0.80, RMSE=150.12kg / hm²) is better than that of the second regressor (R²=0.72, RMSE=179.80kg / hm²), therefore the first regressor is selected.

[0053] The leave-one-year cross-validation strategy is used to find the optimal training set for a single model. When comparing two models, both models undergo leave-one-year cross-validation simultaneously, and evaluation metrics are output to compare the performance of the two models on the same training and validation sets.

[0054] S4: The disaster data is used as the intervention variable, and the time series data and the static data are used as covariates. The causal effect is evaluated based on the intervention variable, covariates and crop yield observations to obtain the crop yield change corresponding to different disaster levels.

[0055] In this embodiment, taking chilling injury as an example, chilling injury is a common type of rice in Northeast China, matching the rice variety described below. The intervention variable is based on the rice low-temperature chilling injury grading standard defined in GB / T 34967-2017, classifying delayed chilling injury into four levels: level 0 indicates no damage, level 1 indicates mild, level 2 indicates moderate, and level 3 indicates severe. In addition, the intervention variable can also be the Palmer Drought Index (PDSI), specifically a calculated index based on a water balance model. Precipitation and temperature data from meteorological station observations are obtained, combined with parameters such as soil available water capacity, and calculated using the standard PDSI calculation formula to obtain continuous monthly PDSI values. Based on the PDSI values, different numerically discrete drought levels are then assigned.

[0056] In one possible approach, a causal forest model is used to assess causal effects, obtaining the changes in crop yield corresponding to different disaster levels, i.e., individual treatment effect values, which also represent the differences in crop yield under different disaster levels compared to no disaster.

[0057] In one possible implementation, the causal forest model is trained using a leave-one-year cross-validation strategy before step S4, during which the validation set corresponding to disaster types that do not appear in the training set is filtered out.

[0058] In this embodiment, a label encoder is used to classify the disaster type of the intervention variable as cold damage. The disaster levels corresponding to cold damage are: no disaster, mild, moderate, and severe. These levels are then converted into numerical values ​​(0, 1, 2, 3). To ensure the reliability of the causal estimation, disaster types that appear in the validation set but not in the training set, specifically cold damage, are filtered out. The time-series data and the static data are used as covariates, and the intervention variable and observed crop yields are input into a causal forest model. The causal forest model fits the influence of one of the intervention variables on the observed crop yield. The influence of the covariates on crop yield is then separated using the causal forest model to obtain the influence of the intervention variable on crop yield; this influence is the individual treatment effect value.

[0059] In one possible implementation, in addition to the individual treatment effect value, the average treatment effect value and the average effect value of the treatment groups can also be obtained. The individual treatment effect value in this application represents the impact of disaster loss on a specific sample; the average treatment effect value represents the impact of the disaster on crop yield loss throughout the entire region, including unaffected individuals; the average effect value of the treatment groups represents the average change in yield caused by the disaster in the sample population that was actually affected by the disaster in the region where the disaster occurred. It is the average of the individual treatment effect values ​​of the intervention in a certain region, i.e., it includes affected individuals.

[0060] In this embodiment, a causal forest model (Causal Forest DML) is used to assess causal effects, with parameters configured as discre_treatment=True, n_estimators=200, min_samples_leaf=15, and random_state=42. The causal forest model is fitted using the training set. For each sample in the validation set, its individual treatment effect value is estimated, ITE_i=Y. 1 _i-Y 0 _i, where Y 0 _i represents the predicted yield under the "no chilling injury" scenario, and Y 1_i represents the predicted yield under the "Level i chilling injury" scenario (i=1, 2, 3). This example focuses on the impact of severe chilling injury (Level 3). The individual treatment effect value ranges from [-650, +50] kg / ha, and the average treatment effect value is the average of the individual treatment effect values ​​of all samples, representing the average impact of chilling injury across the entire region. In this example, the individual treatment effect value is approximately -342 kg / ha, indicating that severe chilling injury results in an average yield reduction of 342 kg / ha. Spatially, the Sanjiang Plain suffers the most severe losses (<-500 kg / ha), while the hilly areas suffer less (≈-200 kg / ha).

[0061] In one possible implementation, S4 also includes a visualization step: plotting the average causal effect curve and the spatial distribution map of individual treatment effects based on the causal effect assessment results.

[0062] In this embodiment, the average yield loss corresponding to different levels of cold damage is displayed by plotting an average causal effect curve, and the spatial heterogeneity of disaster loss is visually displayed by an individual treatment effect spatial distribution map. For example, this embodiment shows that the Sanjiang Plain suffers the most severe loss (<-500 kg / ha), while the hilly areas suffer less (≈-200 kg / ha).

[0063] As a second aspect of this application, a causal inference system for assessing agricultural meteorological disaster losses is also proposed, the system comprising: The data acquisition unit is used to acquire multi-source spatiotemporal data, which includes time-series data, static data, disaster data, and crop yield observations. The time-series data includes meteorological and remote sensing data. As the input end of the system, this unit is responsible for the collection and preliminary structuring of multi-source spatiotemporal data in order to establish a unified spatial-temporal data source and align data from different sources (meteorological stations, remote sensing satellites, soil surveys, and agricultural statistics) to the same sample points and the same time series.

[0064] The preprocessing unit is used to reshape the meteorological and remote sensing data to form a two-dimensional matrix. The two-dimensional matrix is ​​set as the number of samples × the time step. The two-dimensional matrix is ​​then standardized to obtain preprocessed time-series data, which is convenient for the subsequent feature extraction unit to perform time-series processing on the time-series data.

[0065] The feature extraction unit, connected to the data acquisition unit, is used to perform temporal processing on the time-series data to generate deep temporal features. It generates attention weights corresponding to each time step in the time-series data through a multi-head self-attention mechanism. This unit integrates a deep learning model combining a Bidirectional Long Short-Term Memory (BiLSTM) network and a multi-head self-attention mechanism. The BiLSTM is responsible for learning deep dependencies in temporal dynamics from both forward and backward directions, capturing non-linear growth patterns. The self-attention mechanism automatically evaluates and quantifies the importance of each time step (such as different crop growth stages) to the final yield prediction, and its output weights have clear interpretive meaning.

[0066] The dual-path regression unit includes a first regressor and a second regressor, which are connected to the feature extraction unit. The unit is used to input the deep temporal features and the static data into the first regressor for training, and at the same time, to weight the temporal data according to the attention weights and input it into the second regressor along with the static data for training. The interpretation unit uses the SHAP interpreter to interpret and analyze the second regressor, quantifying the impact of environmental characteristics on yield. The causal assessment unit uses disaster data as an intervention variable, time-series data and static data as covariates, and combines them with crop yield observations to assess causal effects and output individual treatment effect values. The individual treatment effect values ​​represent the changes in crop yield corresponding to different disaster levels.

[0067] The dual-path regression unit aims to achieve high-precision crop yield prediction by training two regressors in parallel. Its purpose is to improve prediction accuracy and model interpretability, providing a reliable predictive basis for subsequent causal assessment. The causal forest model, on the other hand, operates independently of the yield prediction process. It directly utilizes observed yield data, intervention variables (disaster severity), and covariates to quantify the causal effects of different disaster severity levels on crop yield, outputting individual treatment effect values ​​to assess the specific impact of disasters on yield. Both models function independently: dual-path regression predicts yield, while causal forest assesses disaster impact, together forming a complete disaster loss assessment process.

[0068] The system achieves accurate assessment of agricultural meteorological disaster losses through a modular architecture. The data acquisition module first collects multi-source spatiotemporal data, including meteorological remote sensing time-series features, static data including soil and topography, discrete numerical values ​​used to evaluate disaster levels, and crop yield observations. The feature extraction module then uses a deep learning network cascaded with LSTM and multi-head self-attention to generate attention weights while extracting time-series features, automatically focusing on key crop growth stages. The dual-path regression module trains two regressors in parallel based on attention weights, using deep time-series features and attention-weighted raw features respectively to predict yield. The optimal model is selected through performance comparison, and the SHAP interpreter is used to interpret and analyze the second regressor, quantifying the impact of environmental features on yield. Finally, the causal assessment module uses a causal forest model to quantify the individual treatment effect value and average treatment effect value of the disaster, generating a spatially heterogeneous visual assessment report, providing a reliable basis for agricultural disaster reduction decision-making.

[0069] As a third aspect of this application, an electronic device is also proposed, including a processor and a memory, the memory storing a computer program, characterized in that the processor executes the computer program to implement the steps of the method described above.

[0070] like Figure 2 As shown, the electronic device is represented in the form of a general-purpose computing device. The components of the electronic device may include, but are not limited to: one or more processors 410, memory 430, and communication bus 440 connecting different system components (including memory 430 and processing unit 410).

[0071] Communication bus 440 represents one or more of several bus architectures, including a memory bus or memory controller, a peripheral bus, a graphics acceleration port, a processor, or a local bus using any of the various bus architectures. For example, these architectures include, but are not limited to, Industry Standard Architecture (ISA) buses, Micro Channel Architecture (MAC) buses, Enhanced ISA buses, Video Electronics Standards Association (VESA) local buses, and Peripheral Component Interconnect (PCI) buses.

[0072] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present disclosure. In this regard, each block in a flowchart or block diagram may represent a unit, segment, or portion of an instruction, which includes one or more executable instructions for implementing a specified logical function. In some alternative implementations, the functions marked in the blocks may occur in a different order than those marked in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, may be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.

[0073] Those skilled in the art will understand that, in the above-described method of the specific implementation, the order in which each step is written does not imply a strict execution order and does not constitute any limitation on the implementation process. The specific execution order of each step should be determined by its function and possible internal logic.

[0074] The various embodiments of this disclosure have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or improvement of the technology in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.

Claims

1. A causal inference method for assessing agricultural meteorological disaster losses, characterized in that, Includes the following steps: S1: Acquire multi-source spatiotemporal data, which includes time-series data, static data, disaster data, and crop yield observations. The time-series data includes meteorological and remote sensing data. S2: Perform time-series processing on the time-series data to generate deep time-series features, and generate attention weights corresponding to each time step in the time-series data through a multi-head self-attention mechanism; S3: Input the deep temporal features and the static data into the first regressor for training. At the same time, weight the temporal data according to the attention weight and input it into the second regressor along with the static data for training. Use the SHAP interpreter to interpret and analyze the second regressor to quantify the impact of environmental features on yield. S4: The disaster data is used as the intervention variable, and the time series data and the static data are used as covariates. The causal effect is evaluated based on the intervention variable, covariates and crop yield observations to obtain the crop yield change corresponding to different disaster levels.

2. The causal inference method for assessing agricultural meteorological disaster losses according to claim 1, characterized in that, The time-step weighting of time-series data based on attention weights described in S3 specifically involves: multiplying the attention weight corresponding to each time step with the observation values ​​of all time-series data within that time step to generate an attention-weighted time-series feature matrix; standardizing the static data to obtain static features, expanding them along the time dimension, and concatenating them with the attention-weighted time-series feature matrix to obtain the attention-weighted raw input; and flattening the attention-weighted raw input along the time-step dimension to form a two-dimensional feature matrix for training the second regressor.

3. The causal inference method for assessing agricultural meteorological disaster losses according to claim 1, characterized in that, The first and second regressors are trained using a one-year cross-validation strategy, and the average of the evaluation metrics from n training iterations is taken as the final evaluation of the model performance; where n is the total number of years.

4. The causal inference method for assessing agricultural meteorological disaster losses according to claim 1, characterized in that, A causal forest model was used to assess causal effects and obtain individual treatment effect values ​​under different disaster levels. These individual treatment effect values ​​represent the changes in crop yield corresponding to different disaster levels.

5. The causal inference method for assessing agricultural meteorological disaster losses according to claim 4, characterized in that, Before step S4, the causal forest model is trained using a one-year leave cross-validation strategy. During the training process, the validation set corresponding to disaster types that do not appear in the training set is filtered out.

6. The causal inference method for assessing agricultural meteorological disaster losses according to claim 1, characterized in that, S4 is followed by a visualization step: drawing the average causal effect curve and the spatial distribution map of individual treatment effects based on the causal effect assessment results.

7. The causal inference method for assessing agricultural meteorological disaster losses according to claim 1, characterized in that, Before S2, the data also includes preprocessing of multi-source spatiotemporal data. The preprocessing includes the following steps: reshaping the meteorological and remote sensing data to form a two-dimensional matrix, the two-dimensional matrix being set as the number of samples × the time step, and standardizing the two-dimensional matrix to obtain preprocessed time series data.

8. The causal inference method for assessing agricultural meteorological disaster losses according to claim 7, characterized in that, Both the first and second regressors are XGBoost regression models.

9. A causal inference system for assessing agricultural meteorological disaster losses, employing a causal inference method for assessing agricultural meteorological disaster losses as described in any one of claims 1 to 8, characterized in that, The system includes: The data acquisition unit is used to acquire multi-source spatiotemporal data, which includes time-series data, static data, disaster data, and crop yield observations. The time-series data includes meteorological and remote sensing data. The feature extraction unit, connected to the data acquisition unit, is used to perform time-series processing on the time-series data to generate deep time-series features, and to generate attention weights corresponding to each time step in the time-series data through a multi-head self-attention mechanism. The dual-path regression unit includes a first regressor and a second regressor, which are connected to the feature extraction unit. The unit is used to input the deep temporal features and the static data into the first regressor for training, and at the same time, to weight the temporal data according to the attention weights and input it into the second regressor along with the static data for training. The interpretation unit uses the SHAP interpreter to interpret and analyze the second regressor, quantifying the impact of environmental characteristics on yield. The causal assessment unit uses disaster data as an intervention variable, time-series data and static data as covariates, and combines them with crop yield observations to assess causal effects and output individual treatment effect values. The individual treatment effect values ​​represent the changes in crop yield corresponding to different disaster levels.

10. A causal inference system for assessing agricultural meteorological disaster losses according to claim 9, characterized in that, include: The preprocessing unit is used to reshape the meteorological and remote sensing data to form a two-dimensional matrix, wherein the two-dimensional matrix is ​​set as the number of samples × the time step, and the two-dimensional matrix is ​​standardized to obtain preprocessed time series data.

11. An electronic device comprising a processor and a memory, the memory storing a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1 to 8.

Images