Rice identification method fusing time sequence spectrum difference and phenological information similarity matching

Abstract

The embodiment of the invention provides a rice identification method fusing time sequence spectrum difference and phenological information similarity matching, and belongs to the technical field of image identification, and the method specifically comprises the steps: combining phenology and land cover products to generate an automatic generation identification sample; training the first sample in each month of the target year in parallel to obtain 12 time sequence extraction classifiers, applying the 12 time sequence extraction classifiers to the second sample to obtain a prediction probability, and dividing the prediction probability into a third sample and a fourth sample; calculating the separable index of each time sequence based on the prediction probability of the third sample, and determining the phenological weight of each phenological period; taking the phenological weight as an adjustment parameter, transferring the prediction probability of the third sample to a time domain, calculating the similarity by using a time-weighted dynamic warping algorithm, cooperatively supervising, learning and training a prediction classifier, and verifying by a fourth sample; and forming an identification model and identifying the target area. Through the scheme disclosed by the invention, the recognition accuracy and adaptability are improved.

Description

Technical Field

[0001] The disclosed embodiments relate to the field of image recognition technology, and in particular to a rice recognition method integrating time-series spectral differences with phenological information similarity matching. Background Art

[0002] At present, rice identification is of great significance for agricultural management and food security. At present, large-scale rice mapping methods using remote sensing technology can be divided into three categories: phenology, machine learning, and the fusion of the two.

[0003] Among them, the phenological method obtains the unique phenological characteristics of the target crop in a single pixel compared with other crops, and identifies rice by adjusting the feature threshold. With the successful launch of high-resolution satellites, time series data has become one of the materials for characterizing the seasonal growth pattern of crops. Some studies have identified rice based on the similarity matching of phenological curves based on a single feature of the time series. However, due to the high requirements of the phenological method on temporal resolution and image availability, rice is difficult to be accurately identified in areas with complex climatic conditions.

[0004] Machine learning algorithms (such as support vector machines, neural networks, etc.) can adapt to more complex features by training samples, which can remove the time restrictions of input images on specific growth stages of crops and are widely used in rice identification. However, the large amount of feature redundancy brought by time series data affects the classification efficiency and accuracy. Traditional feature selection methods ignore the connection between crop stages and the growth laws of the entire phenological period, making it difficult to select important time series spectral features. In addition, feature optimization operations need to be repeated in areas with different rice planting systems.

[0005] To make up for the shortcomings of the above two methods, scholars have integrated the advantages of phenology and machine learning and developed a collaborative algorithm between the two. For example, the existing technology automatically generates samples in small cloudless areas based on phenology, and then transfers machine learning to other cloudy areas. It trains multiple classifiers throughout the season and improves the recognition results of cloudy areas through voting strategies. However, the samples derived based on phenology are affected by the spatiotemporal spectral variability of different geographical environments and observation conditions, resulting in their spatial incompleteness. At the same time, the fine identification of tropical and subtropical rice remains a challenging task. Due to the cloudy and rainy climate in the region, the diverse rice planting systems and complex planting scenarios, it takes a lot of time to zoning according to the planting system. Rice samples have higher requirements on spatial integrity, and the model is difficult to migrate.

[0006] It can be seen that there is an urgent need for a rice identification method with high recognition accuracy and adaptability that integrates time-series spectral differences with phenological information similarity matching. Summary of the invention

[0007] In view of this, the embodiments of the present disclosure provide a rice identification method that integrates time-series spectral differences with phenological information similarity matching, which at least partially solves the problems of poor identification accuracy and adaptability in the prior art.

[0008] The present disclosure provides a rice identification method integrating time series spectral difference and phenological information similarity matching, including:

[0009] Step 1: Obtain phenological parameters based on the ratio model, then combine phenology with rice thematic products to generate high-confidence random rice samples, combine phenology with land cover products to generate random non-rice samples, use morphological filtering operations to purify rice samples and non-rice samples, and divide the purified samples into two non-overlapping parts, namely the first sample and the second sample;

[0010] Step 2: Use the first sample as training data and train it in parallel in each month of the target year to obtain 12 time series extraction classifiers. Apply all classifiers to the second sample to obtain the prediction probability of the time series containing each phenological period, and divide it into the third sample and the fourth sample.

[0011] Step 3, calculating the separability index of each time series based on the predicted probability of the third sample, and determining the phenological weight of each phenological period by normalization;

[0012] Step 4: Using the phenological weight as an adjustment parameter, the prediction probability of the third sample is transferred to the time domain, the similarity is calculated using the time-weighted dynamic warping algorithm, and the prediction classifier is trained by collaborative supervised learning, which is verified by the fourth sample;

[0013] Step 5: The 12 time series extraction classifiers and prediction classifiers are combined into a recognition model to identify the target area and generate rice mapping results.

[0014] According to a specific implementation of the embodiment of the present disclosure, step 1 specifically includes:

[0015] Step 1.1, reconstructing the time series data of rice growth in the sample area, and using the ratio model to analyze the dates corresponding to the key change points in the normalized vegetation index time series curve corresponding to the time series data, and extracting crop phenological parameters, wherein the crop phenological parameters include the start of the growing season, the end of the growing season, the peak of the growing season, the growing season cycle, the kurtosis of the growing season, and the skewness of the growing season;

[0016] Step 1.2, define the preset period before and after the beginning of the growing season as the transplanting period of rice, and perform threshold analysis on the flooding index and the VH band of the Sentinel 1 image during the transplanting period to determine the rice area with flooding signals, while limiting the time period range of the growing season cycle to obtain the potential rice area;

[0017] Step 1.3, using long-term multi-source thematic rice products and phenology-based potential rice for overlay analysis, followed by morphological processing using spatial filtering to obtain a purified potential rice map, and using land cover products to overlay the cultivated land layer to exclude non-cultivated land;

[0018] In step 1.4, a random selection strategy is adopted to generate sample data of rice and other crops and divide them into two non-overlapping parts, namely the first sample and the second sample.

[0019] According to a specific implementation of the embodiment of the present disclosure, step 1.1 specifically includes:

[0020] Step 1.1.1, using the slope extreme value method, in each season, first identify the peak of the normalized vegetation index time series curve as the growing season peak, when the slope is 0 and the time point corresponding to the minimum value of the normalized vegetation index time series curve appears before the growing season peak, it is marked as the beginning of the growing season, and when the slope is 0 and the time point corresponding to the minimum value of the normalized vegetation index time series curve appears after the growing season peak, it is marked as the end of the growing season;

[0021] Step 1.1.2, taking the difference between the beginning of the growing season and the end of the growing season as the growing season period;

[0022] Step 1.1.3, the normalized vegetation index time series curve is regarded as a probability distribution, and the kurtosis of the distribution is regarded as the growing season kurtosis, wherein the expression of the growing season kurtosis is:

[0023]

[0024] Among them, x i Represents the annual daily doy value of the normalized difference vegetation index in the ith month, ndvi i is the normalized difference vegetation index value of the ith month, represents the mean, std represents the standard deviation, and n is a constant of 12;

[0025] Step 1.1.4, the normalized vegetation index time series curve is regarded as a probability distribution, and the skewness of the distribution is regarded as the growing season skewness, wherein the expression of the growing season skewness is:

[0026]

[0027] According to a specific implementation of the embodiment of the present disclosure, step 3 specifically includes:

[0028] The predicted probability sequence in the third sample is subjected to monthly Jeffries-Matus ita distance calculation to obtain a JM time series, which is then normalized to obtain a phenological weight sequence, where the expression of the phenological weight sequence is:

[0029]

[0030] Among them, i represents the i-th time point in the time series, y represents the final phenological weight, and x represents the Jeffries-Matus ita distance;

[0031] The expression of Jeffries-Matusita distance is

[0032] J = 2(1-e -B )

[0033]

[0034] Among them, B represents Bhattacharyya distance, J represents JM distance, m1 and m2 represent the characteristic means of the feature class, and σ1 and σ2 represent the characteristic standard deviation of the feature class.

[0035] According to a specific implementation of the embodiment of the present disclosure, step 4 specifically includes:

[0036] The phenological weight sequence is used as the time weight of the dynamic time planning algorithm. The probability sequences predicted by 12 random forest classifiers on a monthly basis are trained using the K nearest neighbor algorithm in collaboration with the time-weighted dynamic warping algorithm. A prediction classifier for identifying rice is constructed, and the prediction classifier is trained using collaborative supervised learning and verified by the fourth sample.

[0037] The rice identification scheme for integrating time series spectral differences with phenological information similarity matching in the disclosed embodiment includes: step 1, obtaining phenological parameters based on a ratio model, then combining phenology with rice-themed products to generate high-confidence random rice samples, combining phenology with surface cover products to generate random non-rice samples, using morphological filtering operations to purify rice samples and non-rice samples, and dividing the purified samples into two non-overlapping parts, namely, a first sample and a second sample; step 2, using the first sample as training data, training in parallel in each month of the target year, obtaining 12 time series extraction classifiers, and applying all classifiers to the first sample. For the second sample, the predicted probability of the time series including each phenological period is obtained, and it is divided into the third sample and the fourth sample; step 3, based on the predicted probability of the third sample, the separability index of each time series is calculated, and the phenological weight of each phenological period is determined by normalization; step 4, using the phenological weight as an adjustment parameter, the predicted probability of the third sample is transferred to the time domain, the time-weighted dynamic warping algorithm is used to calculate the similarity, and the prediction classifier is trained by collaborative supervised learning, which is verified by the fourth sample; step 5, the 12 time series extraction classifiers and the prediction classifier are formed into a recognition model and the target area is identified to generate the rice mapping result.

[0038] The beneficial effects of the embodiments of the present disclosure are as follows: through the scheme of the present disclosure, the time-weighted dynamic warping algorithm is considered as a bridge, and by coupling phenology and machine learning, the utilization efficiency of time series spectral differences in phenological information matching is improved. The advantages are: automatically constructing high-level features of the corresponding time series from small feature quantities of each phenological period, without the need for decorrelation processing of a large number of features, reducing the dependence on a large number of feature decorrelation operations; using weighted dynamic time warping to achieve feature alignment of phenological information, improving the utilization efficiency of multi-spectral features in phenological information matching. Provide more theoretical solutions for rice mapping in large areas with diverse planting systems and complex scenes. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present disclosure. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0040] Figure 1 A schematic diagram of a flow chart of a rice identification method that integrates time-series spectral differences and phenological information similarity matching provided in an embodiment of the present disclosure;

[0041] Figure 2 A schematic diagram of a framework of a rice identification method that integrates time-series spectral differences and phenological information similarity matching provided by an embodiment of the present disclosure;

[0042] Figure 3 A schematic diagram of a phenological parameter extraction process provided by an embodiment of the present disclosure;

[0043] Figure 4 A schematic diagram of a rice sample generation process provided by an embodiment of the present disclosure;

[0044] Figure 5 A flowchart of a multi-period probability prediction based on machine learning provided in an embodiment of the present disclosure;

[0045] Figure 6 An example schematic diagram of time warping of two rice NDVI curves provided in an embodiment of the present disclosure, wherein (a) represents the time warping matching path using a dynamic programming algorithm; (b) to (c) represent the change from the traditional point pair distance to the dynamic time weighted planning distance;

[0046] Figure 7 A schematic diagram of a remote sensing image data processing flow provided in an embodiment of the present disclosure. DETAILED DESCRIPTION

[0047] The embodiments of the present disclosure are described in detail below with reference to the accompanying drawings.

[0048] The following describes the embodiments of the present disclosure through specific examples, and those skilled in the art can easily understand other advantages and effects of the present disclosure from the contents disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of the present disclosure, rather than all of the embodiments. The present disclosure can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present disclosure. It should be noted that the following embodiments and features in the embodiments can be combined with each other without conflict. Based on the embodiments in the present disclosure, all other embodiments obtained by ordinary technicians in the field without making creative work are within the scope of protection of the present disclosure.

[0049] It should be noted that various aspects of the embodiments within the scope of the appended claims are described below. It should be apparent that the aspects described herein may be embodied in a wide variety of forms, and any specific structure and / or function described herein is merely illustrative. Based on the present disclosure, it should be understood by those skilled in the art that an aspect described herein may be implemented independently of any other aspect, and two or more of these aspects may be combined in various ways. For example, any number of aspects described herein may be used to implement the device and / or practice the method. In addition, other structures and / or functionalities other than one or more of the aspects described herein may be used to implement this device and / or practice this method.

[0050] It should also be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present disclosure. The drawings only show components related to the present disclosure rather than being drawn according to the number, shape and size of components in actual implementation. In actual implementation, the type, quantity and proportion of each component may be changed arbitrarily, and the component layout may also be more complicated.

[0051] Additionally, in the following description, specific details are provided to facilitate a thorough understanding of the examples. However, one skilled in the art will appreciate that the aspects described may be practiced without these specific details.

[0052] The disclosed embodiment provides a rice identification method integrating time series spectral differences with phenological information similarity matching, and the method can be applied to rice mapping in agricultural scenes.

[0053] See also Figure 1 , is a flow chart of a rice identification method that combines time series spectral differences with phenological information similarity matching provided by an embodiment of the present disclosure. Figure 1 and Figure 2 As shown, the method mainly comprises the following steps:

[0054] Step 1: Obtain phenological parameters based on the ratio model, then combine phenology with rice thematic products to generate high-confidence random rice samples, combine phenology with land cover products to generate random non-rice samples, use morphological filtering operations to purify rice samples and non-rice samples, and divide the purified samples into two non-overlapping parts, namely the first sample and the second sample;

[0055] Furthermore, the step 1 specifically includes:

[0056] Step 1.1, reconstructing the time series data of rice growth in the sample area, and using the ratio model to analyze the dates corresponding to the key change points in the normalized vegetation index time series curve corresponding to the time series data, and extracting crop phenological parameters, wherein the crop phenological parameters include the start of the growing season, the end of the growing season, the peak of the growing season, the growing season cycle, the kurtosis of the growing season, and the skewness of the growing season;

[0057] Step 1.2, define the preset period before and after the beginning of the growing season as the transplanting period of rice, and perform threshold analysis on the flooding index and the VH band of the Sentinel 1 image during the transplanting period to determine the rice area with flooding signals, while limiting the time period range of the growing season cycle to obtain the potential rice area;

[0058] Step 1.3, using long-term multi-source thematic rice products and phenology-based potential rice for overlay analysis, followed by morphological processing using spatial filtering to obtain a purified potential rice map, and using land cover products to overlay the cultivated land layer to exclude non-cultivated land;

[0059] In step 1.4, a random selection strategy is adopted to generate sample data of rice and other crops and divide them into two non-overlapping parts, namely the first sample and the second sample.

[0060] Furthermore, the step 1.1 specifically includes:

[0061] Step 1.1.1, using the slope extreme value method, in each season, first identify the peak of the normalized vegetation index time series curve as the growing season peak, when the slope is 0 and the time point corresponding to the minimum value of the normalized vegetation index time series curve appears before the growing season peak, it is marked as the beginning of the growing season, and when the slope is 0 and the time point corresponding to the minimum value of the normalized vegetation index time series curve appears after the growing season peak, it is marked as the end of the growing season;

[0062] Step 1.1.2, taking the difference between the beginning of the growing season and the end of the growing season as the growing season period;

[0063] Step 1.1.3, the normalized vegetation index time series curve is regarded as a probability distribution, and the kurtosis of the distribution is regarded as the growing season kurtosis, wherein the expression of the growing season kurtosis is:

[0064]

[0065] Among them, x i Represents the annual daily doy value of the normalized difference vegetation index in the ith month, ndvi i is the normalized difference vegetation index value of the ith month, represents the mean, std represents the standard deviation, and n is a constant of 12;

[0066] Step 1.1.4, the normalized vegetation index time series curve is regarded as a probability distribution, and the skewness of the distribution is regarded as the growing season skewness, wherein the expression of the growing season skewness is:

[0067]

[0068] In specific implementation, the reconstructed time series data clearly reflects the dynamic process of rice growth and development within a year. By analyzing the dates corresponding to the key change points in the normalized vegetation index time series NDVI curve, a series of crop phenological parameters can be extracted. This paper sets a total of six most critical phenological parameters for extraction. They are the start of the growing season (Start of season, referred to as SOS), the end of the growing season (End of Season, referred to as EOS), the peak of the growing season (Peak of Season, referred to as POS), the cycle of the growing season (Cyc le of Season, referred to as COS), the kurtosis of the growing season (kurtos is of Season, referred to as KOS), and the skewness of the growing season (SOS). For information on phenological parameters, see

[0069] Table 1. Operational procedures Figure 3 shown.

[0070]

[0071] Table 1

[0072] In order to determine SOS and EOS, this paper adopts the slope extreme value method, given that the fixed threshold method may lack flexibility in dealing with changes in different regions and different years, and the seasonal amplitude threshold method is not accurate enough when seasonal changes are drastic. In each season, the peak value of NDVI (i.e., POS) is first identified. When the slope is 0 and the time point corresponding to the minimum NDVI value appears before POS, it is marked as SOS; when the slope is 0 and the time point corresponding to the minimum NDVI value appears after POS, it is marked as EOS. The calculation formula is as follows:

[0073]

[0074] In the formula, T and T2 represent the target time value and the next time value after the target respectively. T Represents the NDVI ratio input at time T. NDVI T is the vegetation index value at time T in the time series. When i = 0, T represents the peak of the growing season. When i = 1, T represents the SOS or EOS of the growing season.

[0075] KOS measures the flatness of the data distribution, that is, a statistic that describes the steepness of the data distribution. For a normal distribution, a data distribution with a large tail has a larger kurtosis value. SOG measures the symmetry of the data, and a value of 0 indicates perfect symmetry. Taking the NDVI time series curve as a probability distribution, the growth process of different crops presents different approximate skewed distributions. Skewness and kurtosis are both descriptions of their phenological information and reflect the statistical characteristics of the rice growth process. The calculation formula is as follows.

[0076]

[0077]

[0078] In the formula, x i Indicates the annual daily doy value of the NDVI curve in the i-th month, ndvi i is the NDVI value of the i-th month, represents the mean, std represents the standard deviation, and n is a constant of 12

[0079] In order to obtain rice samples, the migration method of thematic products requires sufficiently long time series products, and the clustering algorithm has a large randomness. Based on the phenological information obtained from long-term images as the time window, this paper uses the phenological method to automatically generate rice and non-rice crop samples. The samples are purified with the help of existing rice thematic products. Figure 4 The method flow is shown.

[0080] First, the 15 days before and after SOS (a total of one month) are defined as the rice transplanting period. Further, the threshold analysis of the flooding index FLOODS and the VH band of Sentinel 1 during the transplanting period is performed to determine the rice area with flooding signals. At the same time, the COS value is required to be within the range of 100 to 220 days to obtain potential rice areas. The relevant formula is as follows:

[0081] |Ti-ti ndvi |<15; f'(ti ndvi )=0 and f”(ti ndvi )>0 (5)

[0082]

[0083] FLOODS=LSWI Ti(min) -NDVI Ti(min) (7)

[0084] In the formula, ti ndvi It represents the annual daily doy value of the NDVI curve on the time axis, f(ti ndvi ) represents NDVI with respect to ti ndviTi represents the time period of rice transplanting, and VH represents the VH band of Sentinel 1 image. When Rice=1, it means that the area is rice, and when Rice=0, it means that the area is other crops.

[0085] Next, long-term multi-source thematic rice products were used to overlay the potential rice based on phenology. Spatial filtering was then used for morphological processing to obtain a purified potential rice map. In addition, in order to ensure that the negative samples were representative and diverse of other crops and to avoid interference from non-cultivated land on rice mapping, two widely used global land cover products were integrated, including Global Land 30 and ESA World Cover V100, which were used to overlay the cultivated land layer to exclude non-cultivated land.

[0086] Finally, a random selection strategy was adopted to generate sample data of rice and other crops, and divided them into two non-overlapping parts, namely the first sample and the second sample.

[0087] Step 2: Use the first sample as training data and train it in parallel in each month of the target year to obtain 12 time series extraction classifiers. Apply all classifiers to the second sample to obtain the prediction probability of the time series containing each phenological period, and divide it into the third sample and the fourth sample.

[0088] In specific implementation, machine learning methods used for rice mapping include support vector machine (SVM), decision tree (DT), random forest (RF) and neural network (NN). These methods collect images and samples and determine the decision rules for extracting rice based on feature parameters. This paper constructs recognition features from five aspects, including statistical, topographic, spatial, temporal and spectral features. Figure 5 The parameters and method flow of constructing features are shown. The formula of the coefficient of variation feature is as follows:

[0089]

[0090] In the formula, Mean represents the time series mean, and SD represents the standard deviation.

[0091] This paper selects random forest as the classifier for multispectral information extraction. The adjustment of parameters has little effect on the classification results of this method. In addition, the invention is used to evaluate the importance of each input feature in the classification. The key features such as time series and spectrum are analyzed. The first sample is used as training data and trained in parallel in each month of the target year to obtain 12 different classifiers. Finally, these classifiers are applied to the second sample to obtain the predicted probability sequence, which is divided into the third sample and the fourth sample.

[0092] Step 3, calculating the separability index of each time series based on the predicted probability of the third sample, and determining the phenological weight of each phenological period by normalization;

[0093] Based on the above embodiment, step 3 specifically includes:

[0094] The predicted probability sequence in the third sample is subjected to monthly Jeffries-Matus ita distance calculation to obtain a JM time series, which is then normalized to obtain a phenological weight sequence, where the expression of the phenological weight sequence is:

[0095]

[0096] Among them, i represents the i-th time point in the time series, y represents the final phenological weight, and x represents the Jeffries-Matus ita distance;

[0097] The expression of Jeffries-Matusita distance is

[0098] J = 2(1-e -B )

[0099]

[0100] Among them, B represents Bhattacharyya distance, J represents JM distance, m1 and m2 represent the characteristic means of the feature class, and σ1 and σ2 represent the characteristic standard deviation of the feature class.

[0101] During the specific implementation, it was found in the feature importance that different phenological periods had different degrees of influence on the classification results. The credibility of the quantitative prediction probability depends on the ratio of the degree of separability of rice from other crops in the corresponding phenological period. This paper uses the JM distance as the separation index. This method is the most widely recognized and applied separation index calculation method in remote sensing land object classification. It assumes that the probability distribution of land object class features satisfies the normal distribution, and the obtained JM distance ranges from 0 to 2: 0 represents that the two types of land objects cannot be distinguished under this feature, and 2 represents that the two types of land objects can be completely distinguished under this feature condition. In actual situations, there will be no state where JM is equal to 2. It is usually assumed that when JM is greater than 1, the classification work between the corresponding land objects can be carried out. The calculation formula is as follows:

[0102]

[0103] J = 2(1-e -B ) (10)

[0104] In the formula, B represents Bhattacharyya distance, J represents JM distance, m1 and m2 represent the characteristic means of the feature class, and σ1 and σ2 represent the characteristic standard deviation of the feature class.

[0105] This paper calculates the JM distance of the predicted probability series month by month to obtain a JM time series.

[0106] Finally, the normalization process is performed as shown in Formula 11 to obtain the phenological weight sequence.

[0107]

[0108] In the formula, i represents the i-th time point in the time series, y represents the final phenological weight, and x represents the JM distance.

[0109] Step 4: Using the phenological weight as an adjustment parameter, the prediction probability of the third sample is transferred to the time domain, the similarity is calculated using the time-weighted dynamic warping algorithm, and the prediction classifier is trained by collaborative supervised learning, which is verified by the fourth sample;

[0110] Based on the above embodiment, step 4 specifically includes:

[0111] The phenological weight sequence is used as the time weight of the dynamic time planning algorithm. The probability sequences predicted by 12 random forest classifiers on a monthly basis are trained using the K nearest neighbor algorithm in collaboration with the time-weighted dynamic warping algorithm. A prediction classifier for identifying rice is constructed, and the prediction classifier is trained using collaborative supervised learning and verified by the fourth sample.

[0112] In practice, a common task in time series data processing is to compare the similarity of two sequences. However, the lengths of the two time series that need to be compared for similarity may not be equal. In the field of speech recognition, this is manifested in the different speaking speeds of different people. Even if the same person pronounces the same syllable at different times, its duration cannot be exactly the same. In this paper, this feature is manifested in the differences in crop growth cycles. Even the same crop will show different growth rates due to differences in their growing environment. In addition, due to the differences in phenological calendars in different regions, the growth cycle of the same crop is at different time positions within a year. In the extracted high-level feature sequence, this complexity makes it impossible to align features using the Euclidean distance. Therefore, it is necessary to find more flexible methods to evaluate the similarity of time series in order to explain and compare crop growth patterns under different conditions.

[0113] In response to the above problems, dynamic time programming (DTW) describes the time correspondence between the test curve and the reference curve through a time warping function that meets certain conditions, aiming to solve the warping function corresponding to the minimum cumulative distance when the two are matched. The algorithm first uses the distance between each pair of points in the two sequences to calculate the similarity. Then, by performing local scaling on the time axis, the shapes of the two sequences are made as consistent as possible, thereby obtaining the maximum possible similarity. In this article, the predicted probability is a single value, so the Euclidean distance is calculated as the distance between the point pairs. In order to reduce the complexity of the calculation, the dynamic programming DP (dynamic programming) method is used to complete the calculation. Figure 6 The matching path of time warping using dynamic programming algorithm is demonstrated.

[0114] Due to the existence of different planting systems in reality, such as rotation, intercropping and interplanting, the time series data of rice will be disturbed by other crops outside the rice growing period in the same year. For this reason, the phenological weight sequence is used as the time weight of DTW for calculation. The improved method is called TWDTW algorithm. The K nearest neighbor algorithm is used to cooperate with TWDTW for model training to build a classifier for identifying rice.

[0115] Step 5: The 12 time series extraction classifiers and prediction classifiers are combined into a recognition model to identify the target area and generate rice mapping results.

[0116] In specific implementation, 12 time series extraction classifiers and prediction classifiers can be formed into a recognition model to identify the target area and generate rice mapping results. At the same time, the generated rice mapping results can also be returned to the above steps to further train the accuracy and generalization ability of the model.

[0117] The rice identification method provided in this embodiment integrates the time-series spectral difference and the phenological information similarity matching, and considers the time-weighted dynamic warping algorithm as a bridge, and couples phenology with machine learning to improve the utilization efficiency of time-series spectral differences in phenological information matching. The advantages are: automatically constructing high-level features of the corresponding time series from small feature quantities of each phenological period, without the need to decorrelate a large number of features, reducing the dependence on a large number of feature decorrelation operations; using weighted dynamic time warping to achieve feature alignment of phenological information, improving the utilization efficiency of multi-spectral features in phenological information matching. It provides more theoretical solutions for rice mapping in large areas with diverse planting systems and complex scenes.

[0118] The method will be further described below in conjunction with a specific embodiment. In terms of optical data, this paper selects the multispectral image of Sentinel-2. Sentinel-2 consists of two polar-orbiting satellites placed in the same sun-synchronous orbit with a phase of 180° to each other. It uses a multispectral imager (MSI) for surface observation. The 13 spectral bands cover from visible light to near infrared and short-wave infrared.

[0119] For radar data, this paper chooses microwave images provided by Sentinel-1. It is an imaging radar mission that launched satellites in 2014 and can provide continuous all-weather day and night images in the C band. The Sentinel-1 constellation provides high reliability of revisit time, geographic coverage (250 kilometers) and high resolution (10m-40m).

[0120] This paper collects Sentinel-1 and Sentinel-2 images from 2017 to 2023 at intervals of 2 years. Among them, the L2A level product used by Sentinel-2 is obtained by running sen2cor to calculate the L1 level product. Table 2 shows the band information used in the study, and the band with a spatial resolution of 20m is resampled to 10m.

[0121]

[0122] Table 2

[0123] Data quality control is the guarantee for the smooth progress of the experiment. The study area is located in the subtropical region, and the monsoon climate brings abundant precipitation. Cloudy and rainy weather is mostly concentrated in June-October when the southwest monsoon prevails, which is the same period as rice planting and growth, and cloud pollution is inevitable. This paper removes the snow, cloud, and shadow parts of the image through quality control bands, and performs linear interpolation and completion processing on the monthly median composite time series images. Long-term remote sensing images with seamless time and no holes in space at monthly intervals are obtained. This method assumes that there must be cloud-free images at any point in a period of time, but when the cloud cover time is too long, the effect of linear difference noise elimination is not ideal.

[0124] To this end, this paper uses the time domain filtering method Savitzky-Gol ay filter to smooth the sequence. This method has high signal fidelity in areas with stable seasonal cycles such as coniferous forests and farmlands. The SG filter is implemented through the process of convolution, that is, a continuous subset of adjacent data points is fitted with a low-order nonlinear polynomial through the linear least squares method. The advantage is that different window widths can be arbitrarily selected at any position on the same curve to meet the needs of different smoothing filters. In the experiment, the smoothing window in SG is set to 6 and the polynomial order is set to 3.

[0125] The time range of the entire processing process includes the last March of the year before the target year to the first March of the year after the target year. Figure 7 The image data processing flow is shown; Formula 12 is the calculation method for linear difference of time series.

[0126]

[0127] Where t represents the time point when the pixel is missing due to quality issues, and y t is the remote sensing image at time t, and t1 and t2 are the time points of the two known images.

[0128] This paper uses five vegetation indices: NDVI is used to simulate the plant growth process and characterize the vegetation status and canopy structure; EVI compensates for the saturation defect of NDVI under a closed canopy; LSWI is used to enhance the soil moisture content and weaken the information of vegetation and buildings, which has shown good results in transplanting period detection; PSRI reflects the content of vegetation carotenoids and chlorophyll, which can be used to monitor the growth status and aging degree of vegetation; BSI effectively represents the mineral content in the soil. At the same time, a large number of studies have shown that the phenological information provided by remote sensing image time series can greatly improve the classification accuracy and has been widely used in irrigated crop mapping. Therefore, this paper calculates the VI product of the 12-month time series. Table 3 lists the calculation method of the index characteristics.

[0129]

[0130]

[0131] Table 3

[0132] In order to purify the generated samples, they were compared with the rice map obtained by the method disclosed in the present invention. Table 4 summarizes the basic information of Asian rice products used in this paper.

[0133]

[0134] Table 4

[0135] In areas with different land cover types, due to the influence of different objects with the same spectrum and different objects with the same spectrum, it is often easy to introduce false detection when identifying rice. In order to distinguish between rice and other crops, or between rice and natural cover objects, land cover products containing cultivated land information are usually used to remove non-cultivated land to reduce the occurrence of misjudgment. Table 5 lists the global land cover products used in the experiment.

[0136]

[0137]

[0138] Table 5

[0139] Finally, experiments were carried out according to the technical scheme in different cropping system areas in region A. Table 6 lists the accuracy comparison of different rice identification methods.

[0140]

[0141] Table 6

[0142] It can be seen that the PA and UA of the optical-based Eppf-CM are 71.6% and 70.5%, respectively, and OA is 71.1%; DTW based on the time series vegetation index is also a single data source, and the PA and UA have improved and the difference is not large, which are 78.1% and 77.4%, respectively, and OA is 78%. The average PA of the PSPR integrated method based on multi-classifier parallelism has not improved much, and the UA has improved, which are 77.5% and 82%, respectively, and OA is 82%; the average PA and UA of rice in the phenological method Phenology have improved slightly, which are 76.7% and 76.3%, respectively, and OA is 77%; the average PA and UA of rice in the Mf-RF random forest method based on multi-feature selection have further improved, which are 87.6% and 86.4%, respectively, and OA is 86.7%; This shows the importance and development potential of phenological and temporal features in distinguishing rice from other crops. Finally, the average PA and UA of MPSM are 95.2% and 94%, respectively, and OA is 94%. It can be seen that all indicators of MPSM have been greatly improved, and the overall accuracy has increased by nearly 8 percentage points compared with Mf-RF. This shows that MPSM can meet the survey of rice planting distribution at a lower cost. It also shows that coupling time-series spectral differences with crop growth laws can improve the accuracy of rice mapping in large areas with abundant clouds and rain.

[0143] The rice mapping results of this paper are compared with the spatial distribution of existing rice products, showing a high degree of consistency between the two in the spatial distribution of rice. The results of this paper show clearer plots and boundaries, while the existing methods have problems such as unclear plot distinction or failure to separate adjacent plots, and also miss some square areas. It can be seen that the rice identification mapping performed by the method disclosed in this paper has higher accuracy and adaptability.

[0144] It should be understood that various parts of the present disclosure may be implemented in hardware, software, firmware, or a combination thereof.

[0145] The above is only a specific implementation of the present disclosure, but the protection scope of the present disclosure is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present disclosure should be included in the protection scope of the present disclosure. Therefore, the protection scope of the present disclosure should be based on the protection scope of the claims.

Claims

1. A rice identification method integrating time series spectral difference and phenological information similarity matching, characterized in that: include: Step 1: Obtain phenological parameters based on the ratio model, then combine phenology with rice thematic products to generate high-confidence random rice samples, combine phenology with land cover products to generate random non-rice samples, use morphological filtering operations to purify rice samples and non-rice samples, and divide the purified samples into two non-overlapping parts, namely the first sample and the second sample; Step 2: Use the first sample as training data and train it in parallel in each month of the target year to obtain 12 time series extraction classifiers. Apply all classifiers to the second sample to obtain the prediction probability of the time series containing each phenological period, and divide it into the third sample and the fourth sample. Step 3, calculating the separability index of each time series based on the predicted probability of the third sample, and determining the phenological weight of each phenological period by normalization; Step 4: Using the phenological weight as an adjustment parameter, the prediction probability of the third sample is transferred to the time domain, the similarity is calculated using the time-weighted dynamic warping algorithm, and the prediction classifier is trained by collaborative supervised learning, which is verified by the fourth sample; Step 5: The 12 time series extraction classifiers and prediction classifiers are combined into a recognition model to identify the target area and generate rice mapping results.

2. The method according to claim 1, characterized in that , the step 1 specifically includes: Step 1.1, reconstructing the time series data of rice growth in the sample area, and using the ratio model to analyze the dates corresponding to the key change points in the normalized vegetation index time series curve corresponding to the time series data, and extracting crop phenological parameters, wherein the crop phenological parameters include the start of the growing season, the end of the growing season, the peak of the growing season, the growing season cycle, the kurtosis of the growing season, and the skewness of the growing season; Step 1.2, define the preset period before and after the beginning of the growing season as the transplanting period of rice, and perform threshold analysis on the flooding index and the VH band of the Sentinel1 image during the transplanting period to determine the rice area with flooding signals, while limiting the time period range of the growing season cycle to obtain the potential rice area; Step 1.3, using long-term multi-source thematic rice products and phenology-based potential rice for overlay analysis, followed by morphological processing using spatial filtering to obtain a purified potential rice map, and using land cover products to overlay the cultivated land layer to exclude non-cultivated land; In step 1.4, a random selection strategy is adopted to generate sample data of rice and other crops and divide them into two non-overlapping parts, namely the first sample and the second sample.

3. The method according to claim 2, characterized in that , the step 1.1 specifically includes: Step 1.1.1, using the slope extreme value method, in each season, first identify the peak of the normalized vegetation index time series curve as the growing season peak, when the slope is 0 and the time point corresponding to the minimum value of the normalized vegetation index time series curve appears before the growing season peak, it is marked as the beginning of the growing season, and when the slope is 0 and the time point corresponding to the minimum value of the normalized vegetation index time series curve appears after the growing season peak, it is marked as the end of the growing season; Step 1.1.2, taking the difference between the beginning of the growing season and the end of the growing season as the growing season period; Step 1.1.3, the normalized vegetation index time series curve is regarded as a probability distribution, and the kurtosis of the distribution is regarded as the growing season kurtosis, wherein the expression of the growing season kurtosis is: Among them, x i Represents the annual daily doy value of the normalized difference vegetation index in the ith month, ndvi i is the normalized difference vegetation index value of the ith month, represents the mean, std represents the standard deviation, and n is a constant of 12; Step 1.1.4, the normalized vegetation index time series curve is regarded as a probability distribution, and the skewness of the distribution is regarded as the growing season skewness, wherein the expression of the growing season skewness is:

4. The method according to claim 3, characterized in that , the step 3 specifically includes: The predicted probability sequence in the third sample is subjected to monthly Jeffries-Matusita distance calculation to obtain a JM time series, which is then normalized to obtain a phenological weight sequence, where the expression of the phenological weight sequence is: Among them, i represents the i-th time point in the time series, y represents the final phenological weight, and x represents the Jeffries-Matusita distance; The expression of Jeffries-Matusita distance is J=2(1-e -B ) Among them, B represents Bhattacharyya distance, J represents JM distance, m1 and m2 represent the characteristic means of the feature class, and σ1 and σ2 represent the characteristic standard deviation of the feature class.

5. The method according to claim 4, characterized in that , the step 4 specifically includes: The phenological weight sequence is used as the time weight of the dynamic time planning algorithm. The probability sequences predicted by 12 random forest classifiers on a monthly basis are trained using the K nearest neighbor algorithm in collaboration with the time-weighted dynamic warping algorithm. A prediction classifier for identifying rice is constructed, and the prediction classifier is trained using collaborative supervised learning and verified by the fourth sample.

Images