Large-area crop identification method and system based on multi-source space-time distance fusion features

Abstract

The invention provides a large-area crop identification method and system based on multi-source space-time distance fusion features. The method comprises the following steps: firstly, collecting a high-resolution remote sensing image and SRTM elevation data, and constructing a corresponding feature combination variable data set; secondly, aiming at the space-time dynamic change of the growth state of the crops, utilizing the space-time bidirectional characteristic separation index to optimize the characteristics, and capturing the change rules of different crops in the east and west; then, through fusion of collaborative distance features and optimal features, crop growth of different areas is expressed, and feature importance is compared to construct data sets under different fusion features; and finally, integrating the identification probability of each month of the crop growth period by adopting a weaving network month probability random forest classifier, and improving the identification performance by utilizing distance feature interpolation and multi-source feature combination, thereby realizing accurate identification of different crop categories. According to the invention, through feature extraction, optimization and fusion of multi-source data and classifier design, the precision and application effect of large-area crop identification are improved.

Description

Technical Field

[0001] The invention belongs to the technical field of multi-source remote sensing feature fusion and crop identification application, and in particular relates to a large-area crop identification method based on multi-source spatiotemporal distance fusion features. Background Art

[0002] Reliable crop type classification and spatial distribution provide a scientific basis for food security and sustainable agricultural development. The essence of automatic and accurate identification of crop types lies in the spatiotemporal feature expression and complex geographic distance association between different types of crops. At present, the research on automatic sample feature extraction and selection in a large range to achieve accurate mapping of different crops, the development of multi-source medium and high-resolution remote sensing data has promoted the refined monitoring of agriculture in time and space, but ignored the difference in feature distance in the complex crop sowing process. Especially in large-scale regional space, multi-image source feature expression, spatial distance feature fusion and optimized classification algorithm play a key role in determining the accuracy and timeliness of crop classification in large areas. Therefore, the multi-source feature selection and fusion of optical data, radar data and SRTM elevation data face many challenges in crop classification. In crop identification with multi-source remote sensing image feature fusion, multi-source feature fusion provides richer input information for the classifier to improve the accuracy of crop classification. In addition, since different crops have different characteristics and morphologies at different growth stages, and the large differences within the crop class and the small differences between classes are easy to increase the complexity and computational cost of the classifier, it is necessary to enhance the generalization ability of the classification model for crop identification in a large area. Therefore, optimizing feature selection fusion and optimizing classifier design have become urgent problems to be solved in improving crop identification accuracy and automation. To this end, the present invention proposes a large-area crop identification method based on multi-source spatiotemporal distance fusion features. Summary of the invention

[0003] In order to solve the problem of automatic and accurate identification of crops in large areas, the present invention proposes a large-area crop identification method based on multi-source spatiotemporal distance fusion features. The method achieves higher accuracy in identifying different crops and improves the crop identification capability in large areas.

[0004] The above-mentioned object of the present invention can be achieved by the following technical scheme: a large-area crop identification method based on multi-source spatiotemporal distance fusion features, specifically comprising the following steps:

[0005] Step S1, collecting available multi-source image data within the crop growth cycle, extracting feature indexes, and constructing a feature combination variable data set;

[0006] Step S2, feature optimization is performed based on the spatiotemporal bidirectional feature separation index BGSI, and the importance and separation of different crop variables are calculated by the JM distance index. When the JM value is greater than a certain threshold, step S3 is entered;

[0007] Step S3, expressing the preferred spatiotemporal features of different crops within a large area, and combining the distance features DF on the basis of the preferred features to construct an enhanced data set integrating the preferred features and the distance features;

[0008] Step S4, constructing a random forest classifier for the probability of weaving the web month;

[0009] Step S5, input the enhanced data set into the weaving month probability random forest classifier for classification and recognition.

[0010] Furthermore, the multi-source available image data in step S1 include synthetic aperture radar data Sentinel1, optical image data Sentinel2 and SRTM elevation data, and the characteristic index includes optical characteristic index, radar characteristic index and SRTM characteristic index. Then the characteristic combination variable data set includes: the characteristic index combination corresponding to Sentinel2+Sentinel1+SRTM; the characteristic index combination corresponding to Sentinel2+SRTM; the characteristic index combination corresponding to Sentinel1+SRTM.

[0011] Furthermore, the optical characteristic index includes:

[0012] (1) Normalized Difference Vegetation Index:

[0013]

[0014] In the above formula, ρ Nir is the near-infrared band, ρ Red is the red band;

[0015] (2) Renormalized vegetation index:

[0016] REP=705+35*(0.5*(ρ red3 +ρ red )-ρ red1 ) / (ρ red2 -ρ red1 )) / 1000 (2)

[0017] In the above formula, ρ red3 is the red light reflectance, ρ red is the near-infrared light reflectivity, ρ red1 is the reflectivity of the boundary between red light and near-infrared light, which is very sensitive to the chlorophyll concentration of vegetation. red2is the reflectivity of another key point on the red edge, which is used to measure the slope change of the red edge;

[0018] (3) Canopy reflectance index:

[0019]

[0020] In the above formula, ρ Red1 is the red edge light reflectivity, ρ Nir is the near-infrared light reflectance;

[0021] (4) Improved soil adjustment vegetation index 2:

[0022]

[0023] In the above formula, ρ Nir is the near-infrared band, ρ Red is the red band;

[0024] (5) Land surface water index:

[0025]

[0026] In the above formula, ρ Nir is the near-infrared band, ρ swir1 It is the short-wave infrared band;

[0027] (6) Normalized difference tillage index:

[0028]

[0029] In the above formula, ρ swir1 is the near-infrared band, ρ Swir2 It is the red edge band;

[0030] (7) Normalized difference snow index:

[0031]

[0032] In the above formula, ρ Green is the green band, ρ Swir It is the shortwave infrared band.

[0033] Furthermore, radar characteristic indices include:

[0034] (8) Radar polarization characteristics and indices:

[0035] sum=VV+VH (8)

[0036] (9) Radar polarization characteristic difference index:

[0037] subtract=VV-VH (9)

[0038] (10) Radar surface water index:

[0039]

[0040] (11) Radar reflectivity index:

[0041]

[0042] In the above formula, VV means that the radar wave is vertically polarized when transmitting and receiving, that is, it is transmitted vertically upward and received in the vertical direction; VH means that the radar wave is vertically polarized when transmitting and horizontally polarized when receiving, that is, it is transmitted in the vertical direction and received in the horizontal direction.

[0043] Furthermore, the SRTM characteristic index is calculated as follows:

[0044]

[0045] Where aspect is the slope, It is the rate of change of elevation Z in the X direction, that is, the horizontal gradient, and is the elevation difference between adjacent pixels in the X direction in the horizontal distance; It is the rate of change of elevation Z in the Y direction, that is, the vertical gradient, which is the elevation difference between adjacent pixels in the vertical direction; Z is the elevation value of the terrain; X and Y are the horizontal and vertical coordinate axes in the image, expressed in plane coordinates on the ground; atan2 is the inverse tangent function, which is used to calculate the slope angle.

[0046] Furthermore, in step S2, the intra-class and inter-class distances between different crops are first calculated based on the intra-class and inter-class differences, thereby generating a feature separability index feature. The calculation formula of the separability index GSI is as follows:

[0047]

[0048] Among them, Inter-class Distance represents the distance between classes, and Intra-class Distance represents the distance within a class; let the eigenvalue be x t , the category label is y t , time is t, for each category C i , calculate the mean μ of all sample feature values ​​in this category i ;

[0049] Calculate the sum of the squared distances of all samples in the category to the mean, that is, the intra-class variance. The formula is:

[0050]

[0051] The total intra-class distance is the sum of the intra-class distances of all categories:

[0052] Intra-class Distance=∑ i Intra-class Distance i (15)

[0053] Calculate the feature mean μ of all categories i The mean μ of each category is calculated, and the square distance from the mean of each category to the total mean is multiplied by the number of samples in that category. The formula is:

[0054] Inter-class Distance i =|C i |(μ i -μ) 2 (16)

[0055] The total inter-class distance is the sum of the inter-class distances of all categories:

[0056] Inter-class Distance=∑ i Inter-class Distance i (17)

[0057] Then, the GSI is subjected to double row and column iterative normalization using the Sinkhorn-Knopp algorithm. The rows and columns of the matrix are iteratively adjusted and the convergence threshold is set. When the difference between the sum of all rows and columns is close to 1 and is less than the convergence threshold ∈, the algorithm is considered to have converged, and the feature separability index BGSI that can represent different crops is calculated.

[0058] Furthermore, the JM distance is an indicator that measures the separability between different categories in classification problems. The JM distance J between categories a and b is ab It can be expressed as:

[0059]

[0060] Among them, m a and m b are the feature means of categories a and b, respectively, l a and l b is the corresponding characteristic variance, B ab represents the Bhattacharyya distance; J ab Represents the separability measure between sample sets a and b.

[0061] Furthermore, the specific implementation of step S3 is as follows:

[0062] 1) Standardized longitude and latitude data, assuming that the original longitude and latitude data is (lat i ,lon i), the standardization can be done using the following formula:

[0063]

[0064] Among them, μ lat and σ lat are the mean and standard deviation of the latitude data, μ lon and σ lon are the mean and standard deviation of the longitude data respectively;

[0065] 2) Calculate the distance between neighboring points, fixed points or construct a distance matrix as needed; assuming there are n sample points, the Euclidean distance between sample points can be calculated using the following formula:

[0066]

[0067] For each sample point i, find the sample point j that is closest to it:

[0068]

[0069] 3) Use interpolation methods to interpolate unknown locations to obtain more eigenvalues; using the Kriging interpolation method, first you need to establish the covariance function C (h) , where h is the distance between sample points, Kriging interpolation estimates The formula is:

[0070]

[0071] in, is the known eigenvalue of sample point i, λ i is the interpolation weight, obtained by solving the following linear equations:

[0072]

[0073] Here, C(d ij ) is the covariance between sample points i and j, C(d i0 ) is the unknown point x0 and the known point x i The covariance of the distance between them, μ is the Lagrange multiplier;

[0074] 4) Combine the distance feature fusion interpolation result with the original 12 features in step 1 to construct an enhanced feature set; 1i , X 2i , ..., X 12i ] and the standardized longitude and latitude data Distance feature And the interpolation result [K 1i , K 2i , ..., Kni ] to obtain the enhanced distance feature set:

[0075]

[0076] Furthermore, in step S4, the weaving semi-monthly probabilistic random forest classifier is based on the grid, and its input is the spatiotemporal characteristics of the samples in the grid and the training samples from 8 adjacent grid tiles. The grid division is based on the standard grid scheme in a large area, and the study area is divided into 4*4 fixed-size grid units. For each sample point in the grid, the classifier votes on the classification results every several days during the crop growing season and selects the category with the most occurrences as the final category.

[0077] Furthermore, in step S5, the recognition result of the classifier is evaluated and optimized by using a cross-validation technique, and several measurement indicators are calculated in the evaluation and optimization, including using confusion matrix, producer accuracy (PA), user accuracy (UA), overall accuracy (OA), Kappa coefficient and F1-score to evaluate the mapping accuracy;

[0078] Among them, the F1 score calculation formula is defined as follows:

[0079]

[0080] The present invention also provides a large-area crop identification system based on multi-source spatiotemporal distance fusion features, including a processor and a memory, the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the large-area crop identification method based on multi-source spatiotemporal distance fusion features as described in the above technical solution.

[0081] Compared with the prior art, the beneficial technical effects of the present invention are as follows:

[0082] (1) The present invention adopts a data set based on the combination of spectral index features, radar index features and SRTM features. The optimal feature combination in a large area can be determined by constructing a spatiotemporal bidirectional feature separation index (BGSI). The Jeffries-Matusita (JM) distance, an indicator of the separability between different crop features, is used to further determine the importance and separability of the features, thereby ensuring the improvement of the classification accuracy of different crops.

[0083] (2) The present invention further introduces distance features, which are coordinated with preferred features to promote the expression of spatiotemporal features between different types of crops and the association of complex geographical distances of crops in large areas. By coupling the preferred features of multi-source time series images with the spatial distance features of large areas, it is beneficial to the promotion and application of different crop identification in large areas.

[0084] (3) The webbed month probability random forest (WMRF) classification model designed by the present invention performs better than the traditional random forest classifier model in crop classification and identification. Under the condition of complex spatiotemporal heterogeneity of large-scale crops, it verifies the application of multi-source remote sensing features in coordination with spatiotemporal distance features in the identification of different crops in a large area. The fusion of multi-source spatiotemporal distance features has a certain effect on reducing the complexity of the classification model and improving the operating efficiency. The coordination of spatiotemporal distance features highlights the importance of multi-source feature fusion in the classification of complex heterogeneity of large-scale crops, which provides a theoretical basis and support for the classification of crops in large areas.

[0085] In summary, the present invention extracts and optimally combines the features of multi-source remote sensing images, constructs a spatiotemporal bidirectional feature separation index (BGSI), reduces the time for processing feature data, can further determine the importance and separability of optimal features between different crops, and constructs a collaborative spatiotemporal distance feature data set of optimal features to promote the complex geographical distance association of crops in a large area. At the same time, a weaving month probability random forest (WMRF) classification model is designed, which can make full use of the spatiotemporal distance characteristics of multi-source remote sensing images, thereby improving the accuracy and generalization ability of large-scale crop classification and recognition models. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1 is a flow chart of a large-area crop identification method based on multi-source spatiotemporal distance fusion features provided by an embodiment of the present invention;

[0087] Figure 2 This is the distribution map of field sample collection of multiple crops in the Hetao Plain in the study area in 2023;

[0088] Figure 3 It is the average point accuracy and overall accuracy based on the BGSI method integrating multi-source remote sensing feature combination (Sentinel-2, Sentinel-1 and SRTM feature selection);

[0089] Figure 4 It is a comparison of classification accuracy of the Weaving Month Probabilistic Random Forest (WMRF) classification model under different feature fusion datasets;

[0090] Figure 5 It is a large-scale crop mapping under multi-source spatiotemporal distance fusion features and a local crop visualization map corresponding to different feature fusion data sets;

[0091] Figure 6 This is the F1 score result of crop classification under the fusion of multi-source spatiotemporal distance features of different classifiers. DETAILED DESCRIPTION

[0092] In order to facilitate those skilled in the art to understand and implement the present invention, Figure 1It is a flow chart of the technical method. The technical solution of the present invention is further elaborated in detail below in combination with the drawings in the embodiments. It should be understood that the embodiments described here are only used to illustrate and explain a part of the embodiments of the present invention, and are not used to limit the present invention.

[0093] The embodiment of the present invention provides a large-area crop identification method based on multi-source spatiotemporal distance fusion features, which utilizes multi-source remote sensing image feature fusion to improve the timeliness and accuracy of crop identification in a large area. Figure 1 The method provided by the embodiment of the present invention is shown in the flowchart, which specifically includes the following:

[0094] Step 1: Data acquisition and preprocessing: Collect available multi-source image data during the crop growth cycle to extract feature variables and construct a feature combination variable data set;

[0095] The crop sample dataset for field survey in 2023 is divided into 80% of the data as training set and 20% of the data as test set. The time series Sentinel 1 and Sentinel 2 high-resolution remote sensing images and SRTM elevation remote sensing image data of the Hetao Plain from May to October 2023 are preprocessed for different sample areas. Figure 2 This is the distribution map of field sample collection of multiple crops in the Hetao Plain in the study area in 2023;

[0096] Step 2, multi-source data variable acquisition and feature data set construction: multi-source image data (time series optical, time series SAR and SRTM data) are used to extract feature index variables and construct a data set of important combination features. The combined feature data set includes: feature index combination corresponding to Sentinel2+Sentinel1+SRTM; feature index combination corresponding to Sentinel2+SRTM; feature index combination corresponding to Sentinel1+SRTM;

[0097] Step 3: Optimization based on spatiotemporal phenological characteristics of crops: feature extraction of corresponding images of different crops during their growth, optimization based on the spatiotemporal bidirectional feature separation index (BGSI), and on this basis, the importance and separation of different crop variables are calculated using the JM distance index;

[0098] Step 4, construct an enhanced data set under distance feature fusion: Use the preferred feature collaborative distance feature (DF) to construct an enhanced data set under different feature fusions. The feature data set includes: feature index combination corresponding to Sentinel2+Sentinel1+SRTM+spatial distance feature (distance feature fusion DFF); feature index combination corresponding to Sentinel2+Sentinel1+SRTM (multi-source feature fusion FF); feature index combination corresponding to Sentinel2+SRTM; feature index combination corresponding to Sentinel1+SRTM;

[0099] Step 5: Design and optimize the classifier: Use the Web Monthly Probability Random Forest (WMRF) classifier to integrate the recognition probabilities of each month in the crop growth period to achieve recognition of different crops;

[0100] Step 6: The multi-source fusion feature combination dataset is subjected to Weaving Month Probabilistic Random Forest (WMRF) classification and recognition. The recognition performance is verified by comparing the results of different classifiers, and the classification results are generated and the classification recognition accuracy is evaluated to achieve accurate identification of different crops in a large area and crop mapping applications.

[0101] Specifically, the feature acquisition uses multi-source remote sensing images to pre-process the field sample set to obtain sample sets corresponding to different feature indices, and the multi-source feature index is extracted. The specific implementation of step 1 in this embodiment includes the following: The optical feature index is calculated as follows:

[0102] (1) Normalized Difference Vegetation Index:

[0103]

[0104] In the above formula, ρ Nir is the near-infrared band, ρ Red is the red band;

[0105] (2) Renormalized vegetation index:

[0106] REP=705+35*(0.5*(ρ red3 +ρ red )-ρ red1 ) / (ρ red2 -ρ red1 )) / 1000 (2)

[0107] In the above formula, ρ red3 is the red light reflectance, ρ red is the near-infrared light reflectivity, ρ red1 is the reflectivity of the boundary between red light and near-infrared light, which is very sensitive to the chlorophyll concentration of vegetation. red2 is the reflectivity of another key point on the red edge, which is used to measure the slope change of the red edge;

[0108] (3) Canopy reflectance index:

[0109]

[0110] In the above formula, ρ Red1 is the red edge light reflectivity, ρ Nir is the near-infrared light reflectance;

[0111] (4) Improved soil adjustment vegetation index 2:

[0112]

[0113] In the above formula, ρ Nir is the near-infrared band, ρ Red is the red band;

[0114] (5) Land surface water index:

[0115]

[0116] In the above formula, ρ Nir is the near-infrared band, ρ Swir1 It is the short-wave infrared band;

[0117] (6) Normalized difference tillage index:

[0118]

[0119] In the above formula, ρ Swir1 is the near-infrared band, ρ Swir2 It is the red edge band;

[0120] (7) Normalized difference snow index:

[0121]

[0122] In the above formula, ρ Green is the green band, ρ Swir It is the short-wave infrared band;

[0123] The radar signature index is calculated as follows:

[0124] (8) Radar polarization characteristics and indices:

[0125] sum=VV+VH (8)

[0126] In the above formula, VV means that the radar wave is vertically polarized when transmitting and receiving (i.e., transmitting vertically upward and receiving vertically), and VH means that the radar wave is vertically polarized when transmitting and horizontally polarized when receiving (i.e., transmitting vertically and receiving horizontally).

[0127] (9) Radar polarization characteristic difference index:

[0128] subtract=VV-VH (9)

[0129] (10) Radar surface water index:

[0130]

[0131] (11) Radar reflectivity index:

[0132]

[0133] The SRTM characteristic index is calculated as follows:

[0134] (12) Slope:

[0135]

[0136] In the above formula It is the rate of change of elevation (Z) in the X direction (i.e. the gradient in the horizontal direction). It is usually the elevation difference in the X direction between adjacent pixels in the horizontal distance; It is the rate of change of elevation (Z) in the Y direction (i.e., the gradient in the vertical direction), usually the elevation difference between adjacent pixels in the vertical direction in the Y direction; Z is the elevation value of the terrain; X and Y are the horizontal and vertical coordinate axes in the image, usually expressed in plane coordinates on the ground; atan2 is the inverse tangent function, which is used to calculate the slope angle.

[0137] Step 3: The calculation of the bidirectional spatiotemporal feature separation index (BGSI) and the JM distance index is implemented as follows:

[0138] The purpose of feature optimization is to select key features with higher discrimination ability for crop classification, so as to improve classification accuracy and model efficiency. For the bidirectional spatiotemporal feature separation index (BGSI) feature optimization, the intra-class and inter-class distances between different crops are calculated based on the intra-class and inter-class differences, thereby generating a feature separability index (GSI) to quantify the feature's discrimination ability;

[0139] Among them, the calculation formula of feature separability index (GSI) is as follows

[0140]

[0141] Let the eigenvalue be x t , the category label is y t , time is t, for each category C i , i is used to identify each sample and calculate the mean μ of all sample feature values ​​in the category i .

[0142] Calculate the sum of the squared distances of all samples in the category to the mean (intra-class variance), the formula is:

[0143]

[0144] Among them C i represents the sample set of the i-th category. Specifically, t usually represents the index or identifier of a sample, and C i is the set of all samples belonging to category i. So, t refers to the time series data, belonging to category C i a certain moment or a certain sample.

[0145] The total intra-class distance is the sum of the intra-class distances of all categories:

[0146] Intra-class Distance=∑ i Intra-class Distance i (15)

[0147] Calculate the feature mean μ of all categories i The mean μ of each category is calculated, and the square distance from the mean of each category to the total mean is multiplied by the number of samples in that category. The formula is:

[0148] Inter-class Distance i =|C i |(μ i -μ) 2 (16)

[0149] The total inter-class distance is the sum of the inter-class distances of all categories:

[0150] Inter-class Distance=∑ i Inter-class Distance i (17)

[0151] The Sinkhorn-Knopp algorithm is used to perform double row and column iterative normalization on the GSI of different features of specific crops in the same time stage and the GSI of the same feature of specific crops in different time stages to further calculate the importance of spatiotemporal features to the expression of different crops. The Sinkhorn-Knopp algorithm iteratively adjusts the rows and columns of the matrix and sets the convergence threshold. When the difference between the sum of all rows and columns is close to 1 and is less than the convergence threshold ∈, the algorithm is considered to converge. Therefore, GSI is combined with the Sinkhorn-Knopp algorithm to calculate the feature separability index BGSI that can represent different crops. The obtained BGSI matrix is ​​further subjected to feature optimization. According to the spatial dimension of the matrix, features suitable for the specific conditions of the east and west regions can be selected to improve the accuracy and robustness of the classification model. According to the time dimension of the matrix, the key features of different crop growth stages can be determined.

[0152] Jeffries-Matusita (JM) distance is an indicator that measures the separability between different categories in classification problems. Calculating the Jeffries-Matusita (JM) distance can effectively evaluate the separability of different categories in the feature space, thereby solving the problem of quantifying the importance and separability of different crop characteristic variables. For the calculation of the importance and separability of crop characteristic variables in a large area, JM distance is a classic indicator to measure the separability between different categories. It combines the influence of the mean difference and covariance between categories and reflects the degree of overlap of category distribution. The size of the JM distance has the following meanings: the larger the distance, the more separable the categories; the smaller the distance, the more difficult it is to separate the categories. JM distance J of different categories a and b ab It can be expressed as:

[0153]

[0154] Among them, m a and m b are the feature means of categories a and b, respectively, l a and l b is the corresponding characteristic variance, B ab stands for Bhattacharyya distance. ab Represents the separability measure between sample sets a and b. A JM value greater than 1.4 indicates strong separability between different categories. When the strong separability between samples under the feature is verified, an enhanced feature set is constructed.

[0155] Figure 3The mean drop accuracy (MDA) and overall accuracy results are calculated based on the multi-source remote sensing feature combination (Sentinel-2, Sentinel-1 and SRTM feature selection) dataset under the BGSI method. The mean drop accuracy (MDA) is a method used to evaluate the importance of features, especially in random forest models. MDA measures the importance of specific features to the overall model performance. This method helps to identify and select the most influential features, thereby improving the efficiency and accuracy of the model.

[0156] In step 4, the distance matrix of neighboring points and fixed point features is constructed by using the latitude and longitude standardized data; an enhanced feature set is constructed by combining the distance feature fusion Kriging interpolation result with the original 12 features. The specific implementation method is as follows:

[0157] The optimized features and the distance features (DF) will be used to construct an enhanced dataset under the fusion of different features;

[0158] By using the latitude and longitude standardized data to construct the neighboring point distance and fixed point feature distance matrix;

[0159] Standardized longitude and latitude data, assuming the original longitude and latitude data is (lat i ,lon i ), the standardization can be done using the following formula:

[0160]

[0161] Among them, μ lat and σ lat are the mean and standard deviation of the latitude data, μ lon and σ lon are the mean and standard deviation of the longitude data respectively.

[0162] The Kriging interpolation method will be used to interpolate the unknown position to obtain more feature values, and then obtain the distance feature fusion; assuming there are n sample points, the Euclidean distance between the sample points can be calculated using the following formula:

[0163]

[0164] For each sample point i, find the sample point j that is closest to it:

[0165]

[0166] Establish the covariance function C (h) , where h is the distance between sample points, Kriging interpolation estimates The formula is:

[0167]

[0168] Among them, z (xi) is the known eigenvalue of sample point i, λ i To interpolate the weights, we can solve the following linear equations:

[0169]

[0170] Among them, C(d ij ) is the covariance between sample points i and j, C(d i0 ) is the unknown point x0 and the known point x i The covariance of the distance between them, μ is the Lagrange multiplier;

[0171]

[0172] In order to meet the crop growth differences caused by spatial distance in a large area, an enhanced feature set is constructed by merging the distance feature fusion interpolation results with the original 12 features to more comprehensively express the crop characteristics and their regional differences.

[0173] The original 12 features [X 1i , X 2i , ..., X 12i ] and the standardized longitude and latitude data Distance feature And the interpolation result [K 1i , K 2i , ..., K ni ] to obtain the enhanced distance feature set:

[0174]

[0175] The specific implementation of step 5 is as follows:

[0176] A weaving semi-monthly probabilistic random forest classification model (WMRF) was optimized and designed based on grids, with the inputs being the spatiotemporal features of samples within the grid and training samples from eight adjacent grid tiles. The grid division is based on the standard grid scheme within a large area, dividing the study area into 4*4 fixed-size grid cells. For each sample point in the grid, the model votes on the classification results every 15 days during the crop growing season, selects the category with the most occurrences as the final category, and uses the consistency voting mechanism to improve the accuracy of the classification results and conduct accuracy evaluation;

[0177] Among them, the model was optimized and cross-validated, and several accuracy evaluation indicators were calculated, including the use of confusion matrix, producer accuracy (PA), user accuracy (UA), overall accuracy (OA), Kappa coefficient and F1-score to evaluate the mapping accuracy.

[0178] Evaluation of classification mapping accuracy The F1 score calculation formula is defined as follows:

[0179]

[0180] Step 6 is implemented as follows:

[0181] By using different source fusion feature combination data sets, the ability of the Weaving Month Probabilistic Random Forest (WMRF) classifier is verified. By comparing with the results of different classifiers, the crop recognition performance under different multi-source feature fusion data sets is verified, and the accurate recognition of different crops in a large area and crop mapping applications are achieved. In order to evaluate the impact of the WMRF classifier on the classification and recognition capabilities of different crops, it is compared with the accuracy of the RF classification model.

[0182] In order to evaluate the classification and recognition ability of the Weaving Month Probabilistic Random Forest (WMRF) model, the recognition results of different crops under the enhanced dataset containing distance feature fusion constructed in step 4 and the multi-source feature combination dataset are compared. Figure 4 It is a comparison of crop classification accuracy under different multi-source feature fusion datasets of the Weaving Month Probabilistic Random Forest (WMRF) classification model. Figure 5 It is a large-scale crop mapping under multi-source spatiotemporal distance fusion features under the WMRF classifier and a local crop visualization map under different feature fusion data sets.

[0183] Then, we compare the results of different classifiers to verify the recognition performance and compare the classification accuracy of the three classification models:

[0184] (1) Traditional random forest method (RF) model;

[0185] (2) Optimize the partitioned network random forest method (FFRF) model;

[0186] (3) Web-weaving month probabilistic random forest (WMRF) model.

[0187] The three classifiers were compared under the same multi-source spatiotemporal distance feature fusion and optimal parameters. The sample points used in RF are all sample points in the large area, the sample points used in FFRF are the sample points corresponding to each part of the 4*4 block area network, and the sample points of the 3*3 sliding window combined with the adjacent window in the 4*4 block area network used in WMRF. The classification results of crops generated by the three classifiers are compared in Table 1.

[0188] Table 1 Performance comparison of WMRF classification model and other random forest classifier models

[0189]

[0190] In order to further evaluate the classification and recognition accuracy of different classifiers, the effectiveness of the WMRF classifier in large-area crop mapping applications was demonstrated. Figure 6 This is the F1 score result of crop classification under three different classifiers under the fusion of multi-source spatiotemporal distance features.

[0191] An embodiment of the present invention also provides a large-area crop identification system based on multi-source spatiotemporal distance fusion features, including a processor and a memory, the memory being used to store program instructions, and the processor being used to call the stored instructions in the memory to execute the large-area crop identification method based on multi-source spatiotemporal distance fusion features as described in the above technical solution.

[0192] The above-mentioned specific embodiments are merely embodiments of the preferred fusion of corresponding characteristics of several crops of the present invention, and are not intended to limit the present invention. Based on the technical solution of the present invention and the relevant inspiration of the above-mentioned embodiments, equivalent alternative improvements and combinations made by those skilled in the art to the above-mentioned specific embodiments should all be included in the protection scope of the present invention.

Claims

1. A large-area crop recognition method based on multi-source spatiotemporal distance fusion features, characterized in that: The following steps are involved: Step S1, collecting available multi-source image data within the crop growth cycle, extracting feature indexes, and constructing a feature combination variable data set; Step S2, feature optimization is performed based on the spatiotemporal bidirectional feature separation index BGSI, and the importance and separation of different crop variables are calculated by the JM distance index. When the JM value is greater than a certain threshold, step S3 is entered; Step S3, expressing the preferred spatiotemporal features of different crops within a large area, and combining the distance features DF on the basis of the preferred features to construct an enhanced data set integrating the preferred features and the distance features; Step S4, constructing a random forest classifier for the probability of weaving the web month; Step S5, input the enhanced data set into the weaving month probability random forest classifier for classification and recognition.

2. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 1, characterized in that: The multi-source available image data in step S1 include synthetic aperture radar data Sentinel1, optical image data Sentinel2 and SRTM elevation data, and the characteristic index includes optical characteristic index, radar characteristic index and SRTM characteristic index. The characteristic combination variable data set includes: characteristic index combination corresponding to Sentinel2+Sentinel1+SRTM; characteristic index combination corresponding to Sentinel2+SRTM; characteristic index combination corresponding to Sentinel1+SRTM.

3. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 2 is characterized by: Optical characteristic indexes include: (1) Normalized Difference Vegetation Index: In the above formula, ρ Nir is the near-infrared band, ρ Red is the red band; (2) Renormalized vegetation index: REP=705+35*(0.5*(ρ red3 +r red )-r red1 ) / (ρ red2 -r red1 )) / 1000 (2) In the above formula, ρ red3 is the red light reflectance, ρ red is the near-infrared light reflectivity, ρ red1 is the reflectivity of the boundary between red light and near-infrared light, which is very sensitive to the chlorophyll concentration of vegetation. red2 is the reflectivity of another key point on the red edge, which is used to measure the slope change of the red edge; (3) Canopy reflectance index: In the above formula, ρ Red1 is the red edge light reflectivity, ρ Nir is the near-infrared light reflectance; (4) Improved soil adjustment vegetation index 2: In the above formula, ρ Nir is the near-infrared band, ρ Red is the red band; (5) Land surface water index: In the above formula, ρ Nir is the near-infrared band, ρ Swir1 It is the short-wave infrared band; (6) Normalized difference tillage index: In the above formula, ρ Swir1 is the near-infrared band, ρ Swir2 It is the red edge band; (7) Normalized difference snow index: In the above formula, ρ Green is the green band, ρ Swir It is the shortwave infrared band.

4. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 2, characterized in that: Radar signature indices include: (8) Radar polarization characteristics and indices: sum=VV+VH (8) (9) Radar polarization characteristic difference index: subtract=VV-VH (9) (10) Radar surface water index: (11) Radar reflectivity index: In the above formula, VV means that the radar wave is vertically polarized when transmitting and receiving, that is, it is transmitted vertically upward and received in the vertical direction; VH means that the radar wave is vertically polarized when transmitting and horizontally polarized when receiving, that is, it is transmitted in the vertical direction and received in the horizontal direction; The SRTM characteristic index is calculated as follows: Where aspect is the slope, It is the rate of change of elevation Z in the X direction, that is, the horizontal gradient, and is the elevation difference between adjacent pixels in the X direction in the horizontal distance; It is the rate of change of elevation Z in the Y direction, that is, the vertical gradient, which is the elevation difference between adjacent pixels in the vertical direction; Z is the elevation value of the terrain; X and Y are the horizontal and vertical coordinate axes in the image, expressed in plane coordinates on the ground; atan2 is the inverse tangent function, which is used to calculate the slope angle.

5. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 1, characterized in that: In step S2, the intra-class and inter-class distances between different crops are first calculated based on the intra-class and inter-class differences, thereby generating feature separability index features. The calculation formula of the separability index GSI is as follows: Among them, Inter-class Distance represents the distance between classes, and Intra-class Distance represents the distance within a class; let the eigenvalue be x t , the category label is y t , time t, for each category C i , calculate the mean μ of all sample feature values ​​in this category i ; Calculate the sum of the squared distances of all samples in the category to the mean, that is, the intra-class variance. The formula is: The total intra-class distance is the sum of the intra-class distances of all categories: Intra-class Distance=∑ i Intra-class Distance i (15) Calculate the feature mean μ of all categories i The mean μ of each category is calculated, and the square distance from the mean of each category to the total mean is multiplied by the number of samples in that category. The formula is: Inter-class Distance i =|C i |(m i -m) 2 (16) The total inter-class distance is the sum of the inter-class distances of all categories: Inter-class Distance=∑ i Inter-class Distance i (17) Then, the GSI is subjected to double row and column iterative normalization using the Sinkhom-Knopp algorithm. The rows and columns of the matrix are iteratively adjusted and the convergence threshold is set. When the difference between the sum of all rows and columns is close to 1 and is less than the convergence threshold ∈, the algorithm is considered to have converged, and the feature separability index BGSI that can represent different crops is calculated.

6. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 1, characterized in that: The JM distance is an indicator that measures the separability between different categories in a classification problem. The JM distance J between categories a and b ab It can be expressed as: Among them, m a and m b are the feature means of categories a and b, respectively, l a and l b is the corresponding characteristic variance, B ab represents the Bhattacharyya distance; J ab Represents the separability measure between sample sets a and b.

7. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 1, characterized in that: The specific implementation of step S3 is as follows: 1) Standardized longitude and latitude data, assuming that the original longitude and latitude data is (lat i ,lon i ), the standardization can be done using the following formula: Among them, μla t and σ lat are the mean and standard deviation of the latitude data, μ lon and σ lon are the mean and standard deviation of the longitude data respectively; 2) Calculate the distance between neighboring points, fixed points or construct a distance matrix as needed; assuming there are n sample points, the Euclidean distance between sample points can be calculated using the following formula: For each sample point i, find the sample point j that is closest to it: 3) Use interpolation methods to interpolate unknown locations to obtain more eigenvalues; using the Kriging interpolation method, first you need to establish the covariance function C (h) , where h is the distance between sample points, Kriging interpolation estimates The formula is: in, is the known eigenvalue of sample point i, λ i is the interpolation weight, obtained by solving the following linear equations: Here, C(d ij ) is the covariance between sample points i and j, C(d i0 ) is the unknown point x0 and the known point x i The covariance of the distance between them, μ is the Lagrange multiplier; 4) Combine the distance feature fusion interpolation result with the original 12 features in step 1 to construct an enhanced feature set; 1i , X 2i , ..., X 12i ] and the standardized longitude and latitude data Distance feature And the interpolation result [K 1i , K 2i , ..., K ni ] to obtain the enhanced distance feature set:

8. The large-area crop identification method based on multi-source spatiotemporal distance fusion features as claimed in claim 1, characterized in that: In step S4, the weaving semi-monthly probabilistic random forest classifier is based on the grid, and the input is the spatiotemporal characteristics of the samples in the grid and the training samples from 8 adjacent grid tiles. The grid division is based on the standard grid scheme in a large area, and the study area is divided into 4*4 fixed-size grid units. For each sample point in the grid, the classifier votes on the classification results every several days during the crop growing season and selects the category with the most occurrences as the final category.

9. The method for identifying large-area crops based on multi-source spatiotemporal distance fusion features as claimed in claim 1, characterized in that: In step S5, the recognition result of the classifier is evaluated and optimized by using a cross-validation technique. Several measurement indicators are calculated in the evaluation and optimization, including confusion matrix, producer accuracy (PA), user accuracy (UA), overall accuracy (OA), Kappa coefficient and F1-score to evaluate the mapping accuracy. Among them, the F1 score calculation formula is defined as follows:

10. A large-area crop identification system based on multi-source spatiotemporal distance fusion features, characterized by: It includes a processor and a memory, the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the large-area crop identification method based on multi-source spatiotemporal distance fusion features as described in any one of claims 1 to 9.

Images