A quantitative method and system for constructing the influence mechanism of water diversion on NDVI
Through quantitative methods, the problem of incomplete qualitative methods in the existing technology is solved, and the mechanism of impact of water diversion on NDVI is accurately constructed, providing a more accurate reference for engineering practice.
Patent Information
- Application Number
- CN202411568971.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2024-11-05
- Publication Date
- 2025-05-27
- Estimated Expiration
- 2044-11-05
AI Technical Summary
In the prior art, qualitative methods are usually used to analyze the influence mechanism of desert NDVI, and the factors are not comprehensive and the methods are relatively single, so the dual influence mechanism of natural and human factors is still unclear.
A quantitative method is used to construct the mechanism of influence of water diversion on NDVI. By obtaining remote sensing image data from multiple data sources for preprocessing, calculating NDVI values, performing multi-source data fusion, screening the main influence factors, and constructing a set of NDVI influence mechanisms based on linear and nonlinear theories.
A more comprehensive and reliable quantitative analysis of the NDVI impact mechanism is achieved, which eliminates the interference of secondary factors, improves the accuracy of the analysis results, and provides a more accurate reference for engineering practice.
Smart Images

Figure CN119478689B_ABST
Abstract
Description
Technical Field
[0001] The present invention belongs to the technical field of quantitative identification of the attribution of NDVI change, and particularly relates to a quantitative method and system for constructing the influence mechanism of water diversion on NDVI. Background Art
[0002] In recent years, the response mechanism of vegetation growth status to the changing environment has gradually become a research hotspot. Affected by climate change and human activities, the changes in desert vegetation show complexity. Quantitatively attributing the impacts of climate change and human activities on desert vegetation is of great significance for reasonably planning vegetation restoration measures and timely adjusting the scale of ecological water diversion.
[0003] The impacts of desert water diversion on vegetation include natural factors and human factors. Natural factors mainly include rainfall, temperature, wind speed, etc. These factors control the water, heat, and energy conditions for vegetation growth, and thus affect the respiration and photosynthesis of vegetation. Due to the harsh conditions in the Kubuqi Desert, the human factor is mainly the amount of ecological water diversion. NDVI (Normalized Difference Vegetation Index) is an important tool for studying vegetation cover changes and is also the most widely used vegetation index at present, which can better reflect the law of dynamic changes of regional underlying surface vegetation. However, research usually uses qualitative methods to analyze the influence mechanism of desert NDVI, and the factors are not comprehensive and the methods are relatively single. Other quantitative analysis methods for the NDVI influence mechanism have certain limitations and uncertainties, and the dual influence mechanism of natural factors and human factors is not clear and needs to be further improved. How to obtain a more applicable NDVI influence mechanism to provide a more accurate reference basis for engineering practice is particularly important. Summary of the Invention
[0004] In view of this, the present invention provides a quantitative method and system for constructing the influence mechanism of water diversion on NDVI to solve the problems in the prior art that qualitative methods are usually used to analyze the influence mechanism of desert NDVI, and the factors are not comprehensive and the methods are relatively single.
[0005] The technical solution adopted by the present invention is as follows:
[0006] A quantitative method for constructing the influence mechanism of water diversion on NDVI, comprising:
[0007] Step 1: Obtain remote sensing image data of multiple data sources in the region and preprocess each scene of remote sensing image;
[0008] In Step 1, the preprocessing of each scene of remote sensing image specifically includes: geometric correction, radiometric calibration, atmospheric correction, and mosaicking and cropping of each scene of remote sensing image.
[0009] Step 2: Calculate the NDVI value of each remotely sensed image after preprocessing, and calculate the annual average NDVI of each data source in the region based on the NDVI value, respectively obtaining the time series of the average NDVI of each data source in the region;
[0010] In Step 2, the NDVI value of each remotely sensed image is calculated by the following formula:
[0011]
[0012] In the formula, ρ nir is the spectral reflectance of the near-infrared band, and ρ r is the spectral reflectance of the red light band.
[0013] Step 3: (1) Use the interpolation method to fill in the missing average NDVI values in the time series of each data source to process the time series of the average NDVI values of multiple data sources into equal-length sequences; According to Step 2, the multi-year NDVI sequences of different data sources can be obtained. Since the lengths of various NDVI sequences are not equal, and the NDVI change trends retrieved by different satellite remote sensing are consistent, that is, the annual NDVI change curves are similar. Therefore, first use the most complete data (such as MODIS) as the reference data sequence, and the starting year is denoted as t 0 , the data source to be interpolated is Landsat8 (the starting year is denoted as t 1 ) and Sentinel2 (the starting year is denoted as t 2 ), and the common cut-off year of the three types of data is denoted as t n , and MODIS has data at both t 1 and t 2 . The interpolation can be carried out through the following formula.
[0014]
[0015] In the formula, t 0 →t 1 and t 0 →t 2 both represent sequences, and t 1 -t n and t 2 -t n both represent two points. ΔNDVI represents the change amount of NDVI in the corresponding time period.
[0016] (2) Calculate the optimal weight of each data source based on the variance method and spatial characteristics; In (1) above, the different satellite remote sensing data can be processed into equal-length sequences. In order to carry out multi-source satellite data fusion and achieve complementary advantages, the present invention focuses on the data quality and spatial resolution of various satellite data sources, and uses the variance method and assigns weights based on spatial characteristics.
[0017] (2.1) The weights are assigned using the variance method. The weight is the percentage of the variance of each NDVI sequence in the total variance. The product with the smallest variance is given the highest weight, and the product with the largest variance is given the lowest weight. Denote the weight of Sentinel-2 as a 1 , and the weight of Landsat8 as b 1 , and the weight of MODIS as c 1 .
[0018] (2.2) Weights are assigned based on spatial characteristics. Data sources with higher spatial resolution can usually provide more detailed information and can be given higher weights. The weight is the proportion of the spatial resolution of each data source. Denote the weight of Sentinel-2 as a 2 , and the weight of Landsat8 as b 2 , and the weight of MODIS as c 2 .
[0019] Combining the above two methods can consider the main characteristics of various data sources, and the final weights are obtained through arithmetic mean. Denote the weight of Sentinel-2 as a (where ), the weight of Landsat8 as b (where ), and the weight of MODIS as c (where ).
[0020] (3) The corresponding NDVI means within the time series of multiple data sources are fused according to the optimal weights to obtain the NDVI fused data; based on the optimal weights obtained in the above steps (1)-(2), the fused data is calculated according to the following formula.
[0021] NDVI = a·NDVI Sentinel-2 + b·NDVI landsat8 + c·NDVI MODIS (4).
[0022] Step 4: Obtain the alternative influencing factors affecting the NDVI fused data. Take the NDVI fused data as the dependent variable and each alternative influencing factor as the independent variable, and based on correlation test, autocorrelation test, F test, and multicollinearity test, screen out the main influencing factors of the NDVI fused data from the alternative influencing factors;
[0023] The climatic factors affecting desert vegetation are mainly hydrothermal conditions, and the anthropogenic factor is mainly the ecological water diversion volume. Therefore, it is necessary to first determine a sufficient number of factors that can affect vegetation growth. In the present invention, precipitation (Pre), average temperature (Tem), maximum temperature (Tem-max), minimum temperature (Tem-min), wind speed (Win), relative humidity (Rhu), sunshine hours (Ssd), annual water diversion volume (Div), cumulative water diversion volume (Divs), etc. are used as alternative factors. Secondly, taking the NDVI fusion data obtained from Equation (4) as the dependent variable and each alternative factor as the independent variable, multiple-method tests based on the multiple linear regression equation are carried out to gradually screen out the main influencing factors.
[0024] 4.1 The principle of correlation test is as follows:
[0025] The strength of the influence of each influencing factor on the NDVI fusion data is measured by the coefficient of determination of the regression model. The coefficient of determination R 2 The calculation formula is as follows:
[0026]
[0027] In the formula, represents the model prediction value at time i, represents the mean value of the long NDVI sequence, y i represents the measured value of the NDVI fusion data at the i-th moment. Among them, the coefficient of determination R 2 ranges from 0 to 1. The closer the value of the coefficient of determination R 2 is to 1, the stronger the explanatory power of the alternative influencing factor on the NDVI fusion data. The closer the coefficient of determination R 2 is to 0, the weaker the explanatory power of the alternative influencing factor on the NDVI fusion data.
[0028] 4.2 The principle of autocorrelation test is as follows:
[0029] In multiple linear regression, the Durbin-Watson (usually abbreviated as DW) statistic is an index used to test whether there is first-order autocorrelation in the residual sequence. In the present invention, by calculating DW, the autocorrelation of each alternative influencing factor is evaluated. The calculation formula is as follows, and the calculation formula is as follows:
[0030]
[0031] In the formula, e t is the residual of the t-th observation value (i.e., the difference between the measured value and the model prediction value), and n is the total number of observation values;
[0032] Among them, a DW value close to 0 indicates strong positive autocorrelation in the residual sequence, a DW value close to 2 indicates that the residual sequence is close to independent, and a DW value close to 4 indicates strong negative autocorrelation in the residual sequence. Therefore, when DW is close to 2, it indicates that the autocorrelation test is passed.
[0033] 4.3 The principle of the F-test sieve is as follows:
[0034] The present invention uses the F-test for analysis of variance (ANOVA), and the method principle is as follows:
[0035]
[0036] In the formula, SSR is the sum of squared deviations between alternative influencing factors, SSE is the sum of squared deviations within alternative influencing factors, m represents the total number of NDVI fusion data, k represents the number of NDVI fusion data in the training group. Then, the test statistic F is calculated, and the critical value F' corresponding to the significance level of 0.05 is determined. If F > F', the F-test is passed.
[0037] 4.4 The principle of the multicollinearity test is as follows:
[0038] The variance inflation factor is used to measure the collinearity degree of each alternative influencing factor, and the calculation formula is as follows:
[0039]
[0040] In the formula, R 2 is the coefficient of determination. This formula is used to calculate the linear correlation between each alternative factor and other factors, so as to measure the degree of collinearity. When the value of VIF is large, it indicates that there is a strong collinearity between the corresponding independent variable and other independent variables. When the VIF value is less than 10, the multicollinearity test is passed.
[0041] 4.5 Finally, considering the above indicators comprehensively, the main influencing factors that have a greater impact on NDVI are screened out.
[0042] Step 5: Based on the main influencing factors of the screened NDVI fusion data, and combined with linear and nonlinear theories, construct multiple models of the relationship between NDVI and the main influencing factors to form an NDVI influence mechanism set;
[0043] The present invention constructs an NDVI influence mechanism set by combining linear and nonlinear theories. Linear methods include multiple linear regression equations. Nonlinear methods include deep neural networks, random forests, and support vector machines. The principles of various methods are as follows.
[0044] (5.1) Multiple linear regression:
[0045] y = β0 +β 1 X 1 +β 2 X 2 +…+β n X n (9)
[0046] In the formula, X 1 , X 2 , …, X n represent the main influencing factors screened in step 4, β 0 represents the constant of the regression equation, β 1 , β 2 , …, β n represent the optimal weights of the main influencing factors, and y represents the NDVI result predicted by the multiple linear regression model.
[0047] (5.2) Random forest:
[0048] The present invention uses a random forest to construct a non-linear mechanism. Given a training set X = X 1 , X 2 , …, X n and the response variable Y (X is the main influencing factor; Y is NDVI), Bootstrap will repeat B times (b = 1, 2, …, B) of random sampling with replacement to construct the sample set [X b , Y], and then perform tree fitting f b . Before fitting, the training data is projected into a random subspace to increase the variation between different CARTs. To avoid overfitting, the final prediction of RF is defined as the average output of each CART, and the general function is:
[0049]
[0050] In the formula, is the prediction result generated by the CART for the i-th input sample X i in the p-th iteration process.
[0051] (5.3) Support vector machine:
[0052] The present invention uses a support vector machine to construct an influence mechanism. Its basic principle is to find a hyperplane (the optimal hyperplane decision function is as follows), separate various main influencing factors screened in step 4, and maximize the boundary (i.e., the margin) between the two categories, so as to better explain the change of NDVI.
[0053]
[0054] In the formula, Sgn() is the sign function; b* To determine the parameters of the optimal partitioning hyperplane; x, x i ∈R N is an N-dimensional vector, x is a point on the hyperplane, x i is the sample data set, (x·x i ) is the dot product of two vectors; y i ∈{1, 2, ..., k} is the k-class partitioning.
[0055] (5.4) Deep neural network
[0056] The influence mechanism set proposed by the present invention includes a deep neural network. This method mimics the structure and working principle of the human brain neural network. Through hierarchical feature learning and weight adjustment, it can achieve high-performance solutions for complex tasks. Input the preferred influence factors in step 4 into the model to stimulate the operation of the neural network. Each neuron layer receives the output of the previous layer as input and calculates the output through a series of non-linear transformations and weight adjustments. Finally, it is trained by the backpropagation algorithm, that is, by calculating the error between the predicted output and the true output and using the gradient descent method to update the weights and bias values in the network until the network reaches a predetermined performance level, so as to better reveal the changes of NDVI and achieve efficient simulation.
[0057] Step 6: Based on the NDVI fusion data, comprehensively evaluate the performance of various models in the NDVI influence mechanism set by using Pearson correlation coefficient, root mean square error, and mean absolute error indicators to obtain the final NDVI influence mechanism set;
[0058] The present invention divides various models (multiple linear regression, random forest, support vector machine, deep neural network) included in the influence mechanism set into a training set (the first 80% of the data sequence) and a validation set (the last 20% of the data sequence). The performance of the models is evaluated by the Pearson correlation coefficient (PCC), root mean square error (RMSE), and mean absolute error (MAE) in the validation set. The principles of various methods are as follows.
[0059] (6.1) Pearson correlation coefficient (PCC)
[0060] The Pearson correlation coefficient reflects the strength of the linear relationship between the measured value of NDVI and the simulated value of the model. The value range of its absolute value is from 0 to 1. The closer it is to 1, the better the performance of the model.
[0061] Generally, 0.8 < PCC ≤ 1.0 means very strong correlation; 0.6 < PCC ≤ 0.8 means strong correlation; 0.4 < PCC ≤ 0.6 means medium correlation; 0.2 < PCC ≤ 0.4 means weak correlation; 0.0 ≤ PCC ≤ 0.2 means extremely weak or no correlation; PCC ≤ 0.0 means negative correlation.
[0062]
[0063] (6.2) Root Mean Square Error (RMSE)
[0064] The root mean square error is used to evaluate the deviation between the measured NDVI value and the model simulation value. Its value is always non - negative. The smaller the value, the smaller the error, and vice versa, the larger the error.
[0065]
[0066] (6.3) Mean Absolute Error (MAE)
[0067] The mean absolute error (MAE) is used to evaluate the difference between the measured NDVI value and the model simulation value, and measure the magnitude of the average error. The mean absolute error can avoid the problem of error cancellation, and thus can accurately reflect the actual error magnitude.
[0068]
[0069] In formulas (12) - (14), N represents the amount of NDVI data, S i represents the model simulation value at time i, O i represents the measured NDVI value at time i. represents the mean value of the measured NDVI sequence, represents the mean value of the simulated NDVI sequence.
[0070] Step 7: Input the future meteorological data and the corresponding - period NDVI data into the set of NDVI impact mechanisms determined in Step 6, and calculate the future annual water diversion amounts corresponding to different models.
[0071] A quantitative system for constructing the impact mechanism of water diversion on NDVI, comprising:
[0072] A pre - processing module: obtaining remote - sensing image data of multiple data sources in a region, and pre - processing each scene of remote - sensing image;
[0073] A calculation module: calculating the NDVI value of each scene of remote - sensing image after pre - processing, and calculating the annual mean NDVI of each data source in the region based on the NDVI value, respectively obtaining the time series of the mean NDVI of each data source in the region;
[0074] A fusion module: using the interpolation method to fill in the missing mean NDVI values within the time series of each data source, so as to process the time series of the mean NDVI of multiple data sources into an equi - length series;
[0075] Calculating the optimal weight of each data source based on the variance method and spatial characteristics;
[0076] Fuse the corresponding NDVI means within the time series of multiple data sources according to the optimal weights to obtain NDVI fusion data;
[0077] Screening module: Obtain the alternative influencing factors that affect the NDVI fusion data. Take the NDVI fusion data as the dependent variable and each alternative influencing factor as the independent variable, and based on correlation test, autocorrelation test, F test, and multicollinearity test, screen out the main influencing factors of the NDVI fusion data from the alternative influencing factors;
[0078] Model construction module: Based on the main influencing factors of the screened NDVI fusion data, and combined with linear and nonlinear theories, construct models of the relationships between multiple NDVI and the main influencing factors to form a set of NDVI influence mechanisms;
[0079] Evaluation module: Based on the NDVI fusion data, comprehensively evaluate the performance of various models in the set of NDVI influence mechanisms using Pearson correlation coefficient, root mean square error, and mean absolute error indicators to obtain the final set of NDVI influence mechanisms;
[0080] Water diversion volume acquisition module: Input the future meteorological data and the NDVI data for the corresponding period into the set of NDVI influence mechanisms determined by the evaluation module, and calculate the future annual water diversion volumes corresponding to different models.
[0081] In summary, due to the adoption of the above technical solutions, the beneficial effects of the present invention are:
[0082] (1) In the present invention, the advantages of multiple data sources are integrated to provide a more comprehensive and reliable data source for the construction of the NDVI influence mechanism.
[0083] (2) In the present invention, the influencing factors are reasonably screened to remove the interference of secondary factors, thereby ensuring the reliability of the quantitative analysis results.
[0084] (3) The "set of influence mechanisms" constructed by the present invention can combine linear and nonlinear methods to fully explore the relationship between NDVI and its main influencing factors, and has high accuracy.
[0085] (4) In the present invention, the future water diversion range can be obtained according to the set of mechanisms, thereby providing a reference for the development of engineering practices. BRIEF DESCRIPTION OF THE DRAWINGS
[0086] The present invention will be described by way of examples with reference to the accompanying drawings, wherein:
[0087] Figure 1 is the flow chart of the present invention;
[0088] Figure 2 is the flow chart of step 3 of the present invention;
[0089] Figure 3 Schematic diagram of the research area in Embodiment 1 of the present invention;
[0090] Figure 4 Schematic diagram of the interannual change of the NDVI fusion result in Embodiment 1 of the present invention;
[0091] Figure 5 Schematic diagram of the principal component regression test result in Embodiment 1 of the present invention. Detailed implementation manners
[0092] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some but not all of the embodiments of the present invention. Usually, the components of the embodiments of the present invention described and illustrated in the accompanying drawings here can be arranged and designed in various different configurations.
[0093] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed present invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without making creative efforts fall within the scope of protection of the present invention.
[0094] It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other.
[0095] It should be noted that like reference numerals and letters denote like items in the following drawings, and thus, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.
[0096] In the present invention, unless otherwise clearly defined and limited, the first feature being "on" or "under" the second feature may include direct contact between the first and second features, or may include indirect contact between the first and second features through additional features therebetween. Moreover, the first feature being "above", "over" and "on top of" the second feature includes the first feature being directly above and obliquely above the second feature, or merely indicating that the first feature has a higher level height than the second feature. The first feature being "under", "below" and "beneath" the second feature includes the first feature being directly below and obliquely below the second feature, or merely indicating that the first feature has a lower level height than the second feature.
[0097] It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other.
[0098] Embodiment 1
[0099] As Figures 1-4 shown, a quantitative method for constructing the influence mechanism of water diversion on NDVI is disclosed in an embodiment of the present invention, including:
[0100] Step 1: Obtain remote sensing image data of multiple data sources in the region, and preprocess each remote sensing image;
[0101] In Step 1, the preprocessing of each remote sensing image specifically includes: geometric correction, radiometric calibration, atmospheric correction, and mosaicking and cropping of each remote sensing image.
[0102] Step 2: Calculate the NDVI value of each preprocessed remote sensing image, and calculate the annual NDVI mean of each data source in the region based on the NDVI value, respectively obtaining the time series of the NDVI mean of each data source in the region;
[0103] In Step 2, the NDVI value of each remote sensing image is calculated by the following formula:
[0104]
[0105] In the formula, ρ nir is the spectral reflectance of the near-infrared band, and ρ r is the spectral reflectance of the red light band.
[0106] Step 3: (1) Use the interpolation method to fill in the missing NDVI means in the time series of each data source to process the time series of the NDVI means of multiple data sources into an equal-length series; according to Step 2, the multi-year NDVI sequences of different data sources can be obtained. Since the lengths of various NDVI sequences are not equal, and the NDVI change trends retrieved by different satellite remote sensing are consistent, that is, the NDVI interannual change curves are similar. Therefore, first use the most complete data (such as MODIS) as the reference data sequence, and record the starting year as t 0 , the data source to be interpolated is Landsat8 (the starting year is recorded as t 1 ) and Sentinel2 (the starting year is recorded as t 2 ), and the common cut-off year of the three types of data is recorded as t n , and MODIS has data at both t 1 and t 2 . The interpolation can be performed through the following formula.
[0107]
[0108] In the formula, t 0 →t 1 and t 0 →t 2 both represent sequences, t1 -t n and t 2 -t n Both represent two points. ΔNDVI represents the change in NDVI during the corresponding period.
[0109] (2) Calculate the optimal weight for each data source based on the variance method and spatial characteristics; in the above (1), different satellite remote sensing data can be processed into sequences of equal length. To carry out multi-source satellite data fusion and achieve complementary advantages, the present invention focuses on the data quality and spatial resolution of various satellite data sources, and uses the variance method and assigns weights based on spatial characteristics.
[0110] (2.1) Assign weights using the variance method. The weight is the percentage of the variance of each NDVI sequence in the total variance. The product with the smallest variance is given the highest weight, and the product with the largest variance is given the lowest weight. Denote the weight of Sentinel-2 as a 1 , the weight of Landsat8 as b 1 , and the weight of MODIS as c 1 .
[0111] (2.2) Assign weights based on spatial characteristics. Data sources with higher spatial resolution can usually provide more detailed information and can be given higher weights. The weight is the proportion of the spatial resolution of each data source. Denote the weight of Sentinel-2 as a 2 , the weight of Landsat8 as b 2 , and the weight of MODIS as c 2 .
[0112] Combining the above two methods can consider the main characteristics of various data sources, and the final weight is obtained by arithmetic mean. Denote the weight of Sentinel-2 as a (where, ), the weight of Landsat8 as b (where, ), and the weight of MODIS as c (where, ).
[0113] (3) Fuse the corresponding NDVI means within the time series of multiple data sources according to the optimal weight to obtain NDVI fusion data; based on the optimal weight obtained in the above steps (1)-(2), the fusion data is calculated according to the following formula.
[0114] NDVI = a·NDVI Sentinel-2 + b·NDVI landsat8 + c·NDVI MODIS (4).
[0115] Step 4: Obtain the alternative influencing factors affecting the NDVI fusion data. Take the NDVI fusion data as the dependent variable and each alternative influencing factor as the independent variable, and based on correlation test, autocorrelation test, F-test, and multicollinearity test, screen out the main influencing factors of the NDVI fusion data from the alternative influencing factors;
[0116] The main climate factors affecting desert vegetation are hydrothermal conditions, and the main anthropogenic factor is the ecological water diversion volume. Therefore, first determine enough factors that can affect vegetation growth. In the present invention, precipitation (Pre), average temperature (Tem), maximum temperature (Tem-max), minimum temperature (Tem-min), wind speed (Win), relative humidity (Rhu), sunshine hours (Ssd), annual water diversion volume (Div), cumulative water diversion volume (Divs), etc. are used as alternative factors. Secondly, take the NDVI fusion data obtained from Equation (4) as the dependent variable and each alternative factor as the independent variable, and conduct multi-method tests based on the multiple linear regression equation to gradually screen out the main influencing factors.
[0117] 4.1 The principle of the correlation test is as follows:
[0118] Use the coefficient of determination of the regression model to measure the strength of the influence of each influencing factor on the NDVI fusion data. The coefficient of determination R 2 The calculation formula is as follows:
[0119]
[0120] In the formula, represents the model prediction value at the i-th moment, represents the mean value of the long NDVI sequence, and y i represents the measured value of the NDVI fusion data at the i-th moment. Among them, the coefficient of determination R 2 ranges from 0 to 1. The closer the value of the coefficient of determination R 2 is to 1, the stronger the explanatory power of the alternative influencing factor for the NDVI fusion data. The closer the coefficient of determination R 2 is to 0, the weaker the explanatory power of the alternative influencing factor for the NDVI fusion data.
[0121] 4.2 The principle of the autocorrelation test is as follows:
[0122] In multiple linear regression, the Durbin-Watson (usually abbreviated as DW) statistic is an index used to test whether there is first-order autocorrelation in the residual sequence. In the present invention, by calculating DW, the autocorrelation of each alternative influencing factor is evaluated. The calculation formula is as follows, and the calculation formula is as follows:
[0123]
[0124] In the formula, e t is the residual of the t-th observation value (i.e., the difference between the measured value and the model predicted value), and n is the total number of observation values;
[0125] Among them, when the DW value is close to 0, it indicates that there is a strong positive autocorrelation in the residual sequence; when the DW value is close to 2, it indicates that the residual sequence is close to independence; when the DW value is close to 4, it indicates that there is a strong negative autocorrelation in the residual sequence. Therefore, when DW is close to 2, it indicates that the autocorrelation test is passed.
[0126] 4.3 The principle of the F-test screening is as follows:
[0127] The present invention uses the F-test for analysis of variance (ANOVA), and the method principle is as follows:
[0128]
[0129] In the formula, SSR is the sum of squares of deviations between alternative influencing factors, SSE is the sum of squares of deviations within alternative influencing factors, m represents the total number of NDVI fusion data, k represents the number of NDVI fusion data in the training group. Then, the test statistic F is calculated, and the critical value F' corresponding to the significance level of 0.05 is determined. If F > F', then the F-test is passed.
[0130] 4.4 The principle of the multicollinearity test is as follows:
[0131] The variance inflation factor is used to measure the collinearity degree of each alternative influencing factor, and the calculation formula is as follows:
[0132]
[0133] In the formula, R 2 is the coefficient of determination. This formula is used to calculate the linear correlation between each alternative factor and other factors, so as to measure the degree of collinearity. When the value of VIF is large, it indicates that there is a strong collinearity between the corresponding independent variable and other independent variables. When the VIF value is less than 10, the multicollinearity test is passed.
[0134] 4.5 Finally, considering the above indicators comprehensively, the main influencing factors with greater influence on NDVI are screened out.
[0135] Step 5: Based on the main influencing factors of the screened NDVI fusion data, and combining linear and nonlinear theories, construct multiple models of the relationship between NDVI and the main influencing factors to form a set of NDVI influence mechanisms;
[0136] The present invention constructs a set of NDVI impact mechanisms by combining linear and non - linear theories. The linear method includes multiple linear regression equations. The non - linear methods include deep neural networks, random forests, and support vector machines. The principles of each method are as follows.
[0137] (5.1) Multiple linear regression:
[0138] y = β 0 +β 1 X 1 +β 2 X 2 +…+β n X n (9)
[0139] In the formula, X 1 , X 2 , …, X n represent the main influencing factors selected through step 4, β 0 represents the constant of the regression equation, β 1 , β 2 , …, β n represent the optimal weights of each main influencing factor, and y represents the NDVI result predicted by the multiple linear regression model.
[0140] (5.2) Random forest:
[0141] The present invention uses a random forest to construct a non - linear mechanism. Given a training set X = X 1 , X 2 , …, X n and the response variable Y (X is the main influencing factor; Y is NDVI), Bootstrap will repeat B times (b = 1, 2, …, B) of random sampling with replacement to construct the sample set [X b , Y], and then perform tree fitting f b . Before fitting, the training data is projected into a random subspace to increase the variation between different CARTs. To avoid overfitting, the final prediction of RF is defined as the average output of each CART, and the general function is:
[0142]
[0143] In the formula, is the prediction result generated by the CART for the i - th input sample X i during the p - th iteration process.
[0144] (5.3) Support vector machine:
[0145] The present invention uses a support vector machine to construct an impact mechanism. Its basic principle is to find a hyperplane (the optimal hyperplane decision function is as follows), separate various main impact factors screened in step 4, and maximize the boundary (i.e., the margin) between the two categories, so as to better explain the change of NDVI.
[0146]
[0147] In the formula, Sgn() is the sign function; b * is the parameter for determining the optimal partitioning hyperplane; x, x i ∈R N is an N-dimensional vector, x is a point on the hyperplane, x i is the sample data set, (x·x i ) is the dot product of two vectors; y i ∈{1, 2,..., k} is the k-class partitioning.
[0148] (5.4) Deep neural network
[0149] The impact mechanism set proposed by the present invention includes a deep neural network. This method imitates the structure and working principle of the human brain neural network. Through hierarchical feature learning and weight adjustment, it can achieve high-performance solutions for complex tasks. Input the impact factors optimized in step 4 into the model to stimulate the operation of the neural network. Each neuron layer receives the output of the previous layer as input and calculates the output through a series of non-linear transformations and weight adjustments. Finally, it is trained by the backpropagation algorithm, that is, by calculating the error between the predicted output and the true output, and using the gradient descent method to update the weights and bias values in the network until the network reaches the predetermined performance level, so as to better reveal the change of NDVI and thus achieve efficient simulation.
[0150] Step 6: Based on the NDVI fusion data, comprehensively evaluate the performance of various models in the NDVI impact mechanism set using the Pearson correlation coefficient, root mean square error, and mean absolute error indicators to obtain the final NDVI impact mechanism set;
[0151] The present invention divides various models (multiple linear regression, random forest, support vector machine, deep neural network) included in the impact mechanism set into a training set (the first 80% of the data sequence) and a validation set (the last 20% of the data sequence). The performance of the model is evaluated using the Pearson correlation coefficient (PCC), root mean square error (RMSE), and mean absolute error (MAE) in the validation set. The principles of various methods are as follows.
[0152] (6.1) Pearson correlation coefficient (PCC)
[0153] The Pearson correlation coefficient reflects the strength of the linear relationship between the measured NDVI values and the model simulation values. The absolute value of it ranges from 0 to 1. The closer it is to 1, the better the performance of the model.
[0154] Generally, 0.8 < PCC ≤ 1.0 means very strong correlation; 0.6 < PCC ≤ 0.8 means strong correlation; 0.4 < PCC ≤ 0.6 means moderate correlation; 0.2 < PCC ≤ 0.4 means weak correlation; 0.0 ≤ PCC ≤ 0.2 indicates extremely weak or no correlation; PCC ≤ 0.0 indicates negative correlation.
[0155]
[0156] (6.2) Root Mean Square Error (RMSE)
[0157] The root mean square error is used to evaluate the deviation between the measured NDVI values and the model simulation values. Its value is always non - negative. The smaller the value, the smaller the error, and vice versa.
[0158]
[0159] (6.3) Mean Absolute Error (MAE)
[0160] The mean absolute error (MAE) is used to evaluate the difference between the measured NDVI values and the model simulation values, and measure the magnitude of the average error. The mean absolute error can avoid the problem of error cancellation, and thus can accurately reflect the actual error size.
[0161]
[0162] In formulas (12) - (14), N represents the amount of NDVI data, S i represents the model simulation value at time i, and O i represents the measured NDVI value at time i. represents the mean value of the measured NDVI sequence, represents the mean value of the NDVI simulation sequence.
[0163] Step 7: Input the future meteorological data and the corresponding - period NDVI data into the set of NDVI impact mechanisms determined in Step 6, and calculate the future annual water diversion amounts corresponding to different models.
[0164] Example 2
[0165] This example proposes a quantitative system for constructing the impact mechanism of water diversion on NDVI, including:
[0166] Pre - processing module: Obtain the remote - sensing image data of multiple data sources in the region, and pre - process each remote - sensing image;
[0167] Calculation module: Calculate the NDVI value of each remotely sensed image after preprocessing, and calculate the annual average NDVI of each data source in the region based on the NDVI value, respectively obtaining the time series of the average NDVI of each data source in the region;
[0168] Fusion module: Use the interpolation method to fill in the missing average NDVI values in the time series of each data source, so as to process the time series of the average NDVI of multiple data sources into equal-length series;
[0169] Calculate the optimal weight of each data source based on the variance method and spatial characteristics;
[0170] Fuse the corresponding average NDVI values in the time series of multiple data sources according to the optimal weight to obtain the NDVI fusion data;
[0171] Screening module: Obtain the alternative influencing factors affecting the NDVI fusion data, use the NDVI fusion data as the dependent variable, and each alternative influencing factor as the independent variable, and based on the correlation test, autocorrelation test, F test, and multicollinearity test, screen out the main influencing factors of the NDVI fusion data from the alternative influencing factors;
[0172] Model construction module: Based on the main influencing factors of the screened NDVI fusion data, and combined with linear and nonlinear theories, construct models of the relationships between various NDVI and the main influencing factors, forming a set of NDVI influence mechanisms;
[0173] Evaluation module: Based on the NDVI fusion data, comprehensively evaluate the performance of various models in the set of NDVI influence mechanisms using the Pearson correlation coefficient, root mean square error, and mean absolute error indicators to obtain the final set of NDVI influence mechanisms;
[0174] Water diversion volume acquisition module: Input the future meteorological data and the NDVI data in the corresponding period into the set of NDVI influence mechanisms determined by the evaluation module, and calculate the future annual water diversion volume corresponding to different models.
[0175] Example 3
[0176] As Figures 1-4 shown, this example is the specific implementation of Example 1.
[0177] 1. Acquisition and preprocessing of remotely sensed images. In order to respond to the Yellow River ice flood prevention and control and the ecological threat of the Kubuqi Desert, Ordos has explored effective ecological restoration measures since 2014. Part of the ice water and flood water of the Yellow River have been introduced into the low-lying areas on the northern edge of the Kubuqi Desert to relieve the ice flood prevention and control pressure during the ice flood season, and at the same time improve the ecological environment of the northern edge of the Kubuqi Desert (ecological governance area). Based on this background, this example takes the Kubuqi Desert ecological governance area as the research area (the geographical location is as Figure 3As shown). Obtain Sentinel-2, Landsat8, and MODIS satellite images before and after water diversion (from 2004 to 2023), and carry out preprocessing work such as geometric correction, radiometric calibration, atmospheric correction, mosaicking and cropping.
[0178] 2. NDVI data inversion. According to the bands of different remote sensing images, the NDVI is calculated using formula (1) respectively. At the same time, calculate the average NDVI of the vegetation growing season (from June to September) in the study area to obtain the time series of the average NDVI before and after water diversion (from 2004 to 2023). The results are shown in the following table. Among them, the starting time of the Sentinel-2 data (spatial resolution of 10 meters) is 2016, the starting time of the Landsat8 data (spatial resolution of 30 meters) is 2013, and the starting time of the MODIS data (spatial resolution of 1000 meters) is 2004.
[0179] Table 1 Inversion results of multi-source satellite remote sensing NDVI data
[0180]
[0181]
[0182] 3. Multi-source NDVI data fusion. The implementation process of this step is as Figure 2 shown.
[0183] (1) Interpolate the NDVI sequence with equal length. According to step 2, the multi-year NDVI sequences of different data sources can be obtained. Since the time lengths of the data sources are not equal, the NDVI trends and change characteristics inverted by different satellite data sources are basically the same. Therefore, the following method is used for interpolation: First, take the most complete data (such as MODIS) as the reference data sequence, and record the starting year as t 0 , the data source to be interpolated is Landsat8 (the starting year is recorded as t 1 ) and Sentinel2 (the starting year is recorded as t 2 ), and the common cut-off year of the three types of data is recorded as t n , and MODIS has data at both t 1 and t 2 . Interpolation can be carried out through formulas (2)-(3), and the interpolation results are shown in the following table.
[0184] Table 2 Results of equal-length interpolation of NDVI sequence
[0185]
[0186] (2) Calculate the optimal weights using multiple methods. The above steps can process different satellite remote sensing data into sequences of equal length. To carry out multi-source data fusion and achieve complementary advantages, multiple methods are combined to calculate the optimal weights.
[0187] 1) Use the variance method to assign weights. The weight is the percentage of the variance of each NDVI sequence in the total variance. The product with the smallest variance is given the highest weight, and the product with the largest variance is given the lowest weight. Denote the weight of Sentinel-2 as a 1 , the weight of Landsat8 as b 1 , and the weight of MODIS as c 1 .
[0188] 2) Assign weights based on spatial characteristics. Data sources with high spatial resolution can usually provide more detailed information and can be given higher weights. The weight is the proportion of the spatial resolution of each data source. Denote the weight of Sentinel-2 as a 2 , the weight of Landsat8 as b 2 , and the weight of MODIS as c 2 .
[0189] Combining the above two methods can consider the main characteristics of various data sources, and the final weights are obtained through arithmetic mean. Denote the weight of Sentinel-2 as a (where, ), the weight of Landsat8 as b (where, ), and the weight of MODIS as c (where, ). The weight calculation results are shown in the following table.
[0190] Table 3 Calculation results of the optimal weights of each remote sensing data source
[0191]
[0192] (3) Calculate the fused data. Based on the optimal weights obtained from the above steps, the fused data is calculated according to formula (4), and the results are as Figure 4 shown.
[0193] 4. Screen the main influencing factors of NDVI. In this embodiment, the daily-scale meteorological data (precipitation (Pre), average temperature (Tem), maximum temperature (Tem-max), minimum temperature (Tem-min), wind speed (Win), relative humidity (Rhu), sunshine hours (Ssd)) of the Kubuqi Desert Ecological Governance Area from 2004 to 2023, the annual water diversion volume (Div), and the cumulative water diversion volume (Div s)(Desert water diversion factor) and other factors are used as alternative factors. The NDVI fusion data obtained from Equation (3) is used as the dependent variable, and multiple tests (correlation test, autocorrelation test, F-test, multicollinearity test) based on the multiple linear regression equation are carried out using Formulas (4)-(7), so as to gradually screen out the main influencing factors.
[0194] The test results show that when the NDVI fusion data is used as the dependent variable and precipitation, average temperature, wind speed, relative humidity, sunshine hours, and annual water diversion volume are used as independent variables, the correlation coefficient R 2 is 0.87 (relatively high correlation, passing the correlation test); DW is 1.92 (close to 2, passing the autocorrelation test); the significance level of the F-test is lower than 0.01 (passing the F-test); the VIF values of each independent variable are all lower than 10 (passing the multicollinearity test); in addition, the random error term follows a normal distribution (passing the normal distribution test of the random error term and the heteroscedasticity test), and the results are as Figure 5 shown.
[0195] 5. Construct the NDVI influence mechanism. In this embodiment, an NDVI influence mechanism set is constructed based on linear and nonlinear theories (Formulas (8-10)). The linear method includes the multiple linear regression equation. The nonlinear methods include deep neural network, random forest, and support vector machine. Based on the historical data set, the training set and the test set are divided. The first 80% of the data is used as the training set, and the last 20% is used as the test set to continuously optimize the NDVI influence mechanism. The multiple linear regression equation is as follows.
[0196] Y = 0.0126 + 0.3173X 1 - 0.0868X 2 - 0.0698X 3 - 1.0463X 4 + 0.0818X 5 - 0.0406X 6
[0197] In the formula, Y represents NDVI, and X 1 represents annual precipitation, X 2 represents average temperature, X 3 represents wind speed, X 4 represents relative humidity, X 5 represents sunshine hours, X 6 represents annual water diversion volume.
[0198] Finally, Formulas (12)-(14) are used to evaluate the effects of various models in the influence mechanism set, and the results are shown in Table 4. Among them, if PCC ≥ 0.6, RMSE ≤ 0.2, and MAE ≤ 0.2 are satisfied simultaneously, it means that the simulation effect of the model is good and can be further used for the prediction of water diversion volume.
[0199] Table 4 Evaluation Results of NDVI Impact Mechanism Set
[0200]
[0201] Based on the calculation results in Table 4 and combined with the above evaluation criteria, deep neural network and random forest are selected for subsequent water diversion volume prediction.
[0202] 6. Adjust the future water diversion area. Research shows that when 0 < NDVI ≤ 0.2, it belongs to the low vegetation coverage area; when 0.2 < NDVI ≤ 0.4, it belongs to the medium-low coverage area; when 0.4 < NDVI ≤ 0.6, it belongs to the medium vegetation coverage area; when 0.6 < NDVI ≤ 0.8, it belongs to the medium-high vegetation coverage area; when 0.8 < NDVI ≤ 1, it belongs to the high coverage area. In this embodiment, considering the climate conditions of the desert and the growth characteristics of vegetation comprehensively, taking the local vegetation reaching the medium coverage level (NDVI: 0.4 - 0.6) as the future planning goal, and at the same time inputting the prediction results of PRECIS (Providing Regional Climates for Impacts Studies) for each meteorological element (precipitation, average temperature, relative humidity, sunshine hours, wind speed) from 2025 to 2034 into the mechanism set determined in step 6, calculating the annual water diversion volume and comprehensive water diversion situation corresponding to various models from 2025 to 2050. The results are shown in Table 5, and further set the future water diversion volume range to be 31.78 - 78.53 million m 3 , and the results can provide reference for engineering practice.
[0203] Table 5 Optimal Water Diversion Volume Prediction Results from 2025 to 2050 (Water Diversion Volume (10,000 m 3 ))
[0204]
[0205]
[0206] The circuits, electronic components and modules involved are all prior arts, which can be fully realized by those skilled in the art without further elaboration. The content protected by the present invention does not involve the improvement of software and methods either.
[0207] The various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. For the same or similar parts among the various embodiments, reference can be made to each other.
[0208] The foregoing description of the disclosed embodiments enables those skilled in the art to practice or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.
Claims
1. A quantitative method for constructing the mechanism of water diversion's impact on NDVI, characterized in that: include: Step 1: Obtain remote sensing image data from multiple data sources in the region and pre-process each remote sensing image; Step 2: Calculate the NDVI value of each remote sensing image after preprocessing, and calculate the annual NDVI mean of each data source in the region based on the NDVI value, and obtain the time series of the NDVI mean of each data source in the region; Step 3: Use the interpolation method to fill in the missing NDVI mean values in the time series of each data source, so as to process the time series of NDVI mean values of multiple data sources into sequences of equal length; Calculate the optimal weight of each data source based on variance method and spatial characteristics; According to the optimal weight, the corresponding NDVI means in the time series of multiple data sources are fused to obtain NDVI fused data; Step 4: Obtain the alternative influencing factors that affect the NDVI fusion data, take the NDVI fusion data as the dependent variable, and each alternative influencing factor as the independent variable, and screen out the main influencing factors of the NDVI fusion data from the alternative influencing factors based on correlation test, autocorrelation test, F test and multicollinearity test; Step 5: Based on the main influencing factors of the screened NDVI fusion data, a variety of models of the relationship between NDVI and the main influencing factors are constructed in combination with linear and nonlinear theories to form a set of NDVI influencing mechanisms; Step 6: Based on the NDVI fusion data, the Pearson correlation coefficient, root mean square error, and mean absolute error indicators are used to comprehensively evaluate the performance of various models of the NDVI impact mechanism set to obtain the final NDVI impact mechanism set; Step 7: Input the future meteorological data and the NDVI data of the corresponding period into the NDVI influencing mechanism set determined in step 6, calculate the future annual water diversion corresponding to different models, and comprehensively consider the results of multiple models to derive the water diversion range.
2. A quantitative method for constructing the mechanism of water diversion's impact on NDVI according to claim 1, characterized in that: In step 1, the preprocessing of each remote sensing image specifically includes: Each remote sensing image is subjected to geometric correction, radiation calibration, atmospheric correction, and stitching and cropping.
3. A quantitative method for constructing the mechanism of water diversion's impact on NDVI according to claim 1, characterized in that: In step 2, the NDVI value of each remote sensing image is calculated by the following formula: In the formula, ρ nir is the spectral reflectance in the near-infrared band, ρ r is the spectral reflectance in the red light band.
4. A quantitative method for constructing the mechanism of water diversion's impact on NDVI according to claim 1, characterized in that: In step 4, the main influencing factors of NDVI fusion data were screened out based on correlation test, including: The determination coefficient of the regression model is used to measure the influence of each influencing factor on the NDVI fusion data. 2 The calculation formula is as follows: In the formula, represents the model prediction value at time i, The mean of the long series of NDVI, y i represents the measured value of NDVI fusion data at the i-th moment, where the determination coefficient R 2 The value range is between 0 and 1, and the determination coefficient R 2 The closer the value is to 1, the stronger the explanatory power of the alternative influencing factors on the NDVI fusion data. 2 The closer it is to 0, the weaker the explanatory power of the alternative influencing factors on the NDVI fusion data.
5. A quantitative method for constructing the mechanism of water diversion's impact on NDVI according to claim 1, characterized in that: In step 4, the main influencing factors of NDVI fusion data were screened out based on the autocorrelation test, including: By calculating DW, the autocorrelation of each candidate influencing factor is evaluated. The calculation formula is as follows: In the formula, e t is the residual of the t-th observation, and n is the total number of observations; Among them, a DW value close to 0 indicates that the residual sequence has a strong positive autocorrelation, and a DW value close to 2 indicates that the residual sequence is close to independence.
6. A quantitative method for constructing the mechanism of water diversion's impact on NDVI according to claim 1, characterized in that: In step 4, the main influencing factors of NDVI fusion data were screened out based on the F test, including: The F test was used for variance analysis as follows: In the formula, SSR is the sum of squares of deviations between the alternative influencing factors, SSE is the sum of squares of deviations within the alternative influencing factors, m represents the total number of NDVI fusion data, and k represents the number of NDVI fusion data in the training group. Then, the test statistic F is calculated, and the critical value F' corresponding to the significance level of 0.05 is determined. If F>F', the F test is passed.
7. A quantitative method for constructing the mechanism of the impact of water diversion on NDVI according to claim 1, characterized in that: In step 4, the main influencing factors of NDVI fusion data were screened out based on multicollinearity test, including: The variance inflation coefficient is used to measure the collinearity of each candidate influencing factor. The calculation formula is as follows: In the formula, R 2 It is the determination coefficient. When the VIF value is large, it indicates that there is strong collinearity between the corresponding independent variable and other independent variables. When the VIF value is less than 10, the multicollinearity test is passed.
8. A quantitative method for constructing the mechanism of water diversion's impact on NDVI according to claim 1, characterized in that: In step 5, the models of the relationship between the various NDVIs and the main influencing factors specifically include the following models: multiple linear regression equation, deep neural network, random forest and support vector machine.
9. A quantitative system for constructing the mechanism of water diversion's influence on NDVI, used to implement the quantitative method for constructing the mechanism of water diversion's influence on NDVI as claimed in claim 1, characterized in that: include: Preprocessing module: obtain remote sensing image data from multiple data sources in the region and preprocess each remote sensing image; Calculation module: Calculate the NDVI value of each remote sensing image after preprocessing, and calculate the annual NDVI mean of each data source in the region based on the NDVI value, and obtain the time series of the NDVI mean of each data source in the region; Fusion module: Use interpolation method to fill in the missing NDVI mean values in the time series of each data source, so as to process the time series of NDVI mean values of multiple data sources into equal-length sequences; Calculate the optimal weight of each data source based on variance method and spatial characteristics; According to the optimal weight, the corresponding NDVI means in the time series of multiple data sources are fused to obtain NDVI fused data; Screening module: Obtain the alternative influencing factors that affect the NDVI fusion data, take the NDVI fusion data as the dependent variable, and each alternative influencing factor as the independent variable, and screen out the main influencing factors of the NDVI fusion data from the alternative influencing factors based on correlation test, autocorrelation test, F test and multicollinearity test; Model building module: Based on the main influencing factors of the screened NDVI fusion data, a variety of models of the relationship between NDVI and the main influencing factors are constructed in combination with linear and nonlinear theories to form a set of NDVI influencing mechanisms; Evaluation module: Based on NDVI fusion data, the Pearson correlation coefficient, root mean square error, and mean absolute error indicators are used to comprehensively evaluate the performance of various models of the NDVI impact mechanism set to obtain the final NDVI impact mechanism set; Water diversion acquisition module: input future meteorological data and NDVI data of the corresponding period into the NDVI impact mechanism set determined by the evaluation module, and calculate the future annual water diversion corresponding to different models.
Citation Information
Patent Citations
A method for predicting concrete durability based on data mining and artificial intelligence algorithm
AU2020101854A4
Vegetation change cause recognition method considering spatial correlation
CN112907113A