Temporary land interpretation and reclamation verification method and system

Abstract

The invention belongs to the technical field of natural resource management, and discloses a temporary land interpretation and reclamation verification method and system, and the method comprises the steps: employing a GEE cloud computing platform and a Sentinel-2 satellite remote sensing image, and calculating and generating a typical time sequence characteristic curve of the remote sensing characteristic indexes of a temporary land and a non-temporary land; calculating a time sequence characteristic curve of the remote sensing characteristic indexes of the land parcel to be judged, calculating DTW similarity distances between the time sequence characteristic curve of the land parcel to be judged and typical time sequence characteristic curves of the temporary land and the non-temporary land, and judging whether the land parcel is the temporary land; and for the temporary land at the reclamation period, calculating the slope of the time sequence characteristic curve according to the remote sensing characteristic indexes at the use period and the reclamation period, and judging whether the land reclamation of the temporary land is completed or not. According to the method, the time cost and the labor cost for manually and visually interpreting and checking the land reclamation condition after the temporary land use period is up can be greatly reduced, and the government supervision level and the decision analysis capability are improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of natural resource management, and relates to a temporary land use interpretation and reclamation verification method and system, and in particular to a technology and method for monitoring and supervising temporary land use based on satellite remote sensing data. Background Art

[0002] Temporary land is an important form of land use in my country. It is state-owned land or land collectively owned by farmers that needs to be temporarily used for construction projects and geological surveys. It provides great convenience for engineering construction and project construction, and plays a positive role in ensuring construction projects and geological surveys. The use period of temporary land is generally within 2 years, and it can be 4 years in special circumstances. For temporary land that occupies cultivated land, forest land, and agricultural land, if the land is not reclaimed in time after the expiration of the use period, it will cause land resource loss, cultivated land damage and other problems. The traditional temporary land supervision method mainly relies on land inspection and law enforcement. After on-site investigation, law enforcement opinions are issued to the project user or the competent department, and then the competent department supervises and rectify. This method makes it difficult to comprehensively, timely, accurately and objectively supervise temporary land projects.

[0003] Satellite remote sensing has the characteristics of fast response, wide range, short cycle, strong timeliness, and repeatable observation. It can continuously track, monitor over a large area, and accurately identify various types of ground objects, providing data support for the monitoring and supervision of temporary land use. However, due to the large differences in the use cycle of temporary land and the complexity of its types, the existing use of remote sensing technology for temporary land use interpretation and reclamation verification is generally judged through manual visual interpretation, and the automation of temporary land use interpretation and land reclamation verification has not yet been achieved.

[0004] Through the above analysis, the problems and defects of the existing technology are as follows: Due to the large differences in the reclamation cycle of temporary land and the complexity of types, the existing remote sensing temporary land interpretation and reclamation verification are judged by manual visual interpretation and have not yet been automated. The traditional temporary land supervision method is difficult to supervise temporary land projects comprehensively, timely, accurately and objectively. The manual visual remote sensing interpretation method requires a lot of manpower and time costs, and is easily affected by human subjective factors, resulting in inconsistent standards and errors. Summary of the invention

[0005] In order to overcome the difficulty of land reclamation verification in temporary land supervision, the disclosed embodiment of the present invention provides a temporary land interpretation and reclamation verification method and system, and in particular, relates to a temporary land interpretation and reclamation verification method and system based on time-series satellite remote sensing images. The technical solution is as follows:

[0006] The present invention is implemented in this way: a temporary land use interpretation and reclamation verification method comprises:

[0007] S1. Using the GEE cloud computing platform and Sentinel-2 satellite remote sensing images, calculate and generate sample characteristics of typical time series characteristic curves of remote sensing characteristic indexes of temporary land and non-temporary land;

[0008] S2, calculating the time series characteristic curve of the remote sensing characteristic index of the plot to be judged, and after linear interpolation or linear compression, calculating the DTW similarity distance between the time series characteristic curve of the remote sensing characteristic index of the plot to be judged and the typical time series characteristic curves of temporary land and non-temporary land, and judging whether the plot is temporary land based on the DTW distance criterion;

[0009] S3. For temporary land that has expired, the slope of the time series characteristic curve is calculated based on the remote sensing characteristic index at the time of expiration of the use period and the time of expiration of the reclamation period to determine whether the temporary land has completed land reclamation.

[0010] Before step S1, the Sentinel-2 satellite remote sensing images are pre-processed by image screening, cloud removal and image cropping; temporary land and non-temporary land samples are obtained to generate time series information of the entire life cycle of temporary land; Sentinel-2 satellite remote sensing image screening is based on the GEE platform to set the date according to the temporary land vector range and the approval time and expiration time, and the image data with cloud cover less than 20% during the temporary land use period and the land reclamation period are screened, and the obtained Sentinel-2 satellite remote sensing images are de-clouded using the QA60 band.

[0011] In step S1, sample features of typical time series characteristic curves of remote sensing characteristic indexes of temporary land and non-temporary land are calculated and generated, including:

[0012] According to the vector range shp file and kml file of temporary land, the GEE cloud platform and Sentinel-2 satellite remote sensing images are used to obtain the time series NDVI, NDBI and NDSI remote sensing characteristic indexes of temporary land during the land use period and the full year time series NDVI, NDBI and NDSI remote sensing characteristic indexes of non-temporary land.

[0013] By performing maximum and minimum normalization processing on the time series characteristic curve, the time series MMNDVI, MMNDBI and MMNDSI values ​​of typical temporary land and non-temporary land patches were obtained;

[0014] The sample characteristics of typical time series characteristic curves of three remote sensing characteristic indices, including MMNDVI, MMNDBI and MMNDSI, for temporary land use and non-temporary land use categories were obtained by mean calculation.

[0015] Furthermore, the time series NDVI, NDBI and NDSI remote sensing characteristic indexes of temporary land use during the land use period are obtained, including:

[0016] The Sentinel-2 satellite sensor includes 12 spectral bands, and the normalized vegetation index NDVI, normalized building index NDBI and normalized soil index NDSI are calculated by band operation to obtain the NDVI, NDBI and NDSI values ​​of satellite remote sensing data at different phases.

[0017] Among them, the calculation expressions of NDVI, NDBI and NDSI are as follows:

[0018]

[0019] Where NIR is the surface reflectance of the near infrared band of the Sentinel-2 satellite remote sensing image, RED is the surface reflectance of the red band of the Sentinel-2 satellite remote sensing image, SWIR is the surface reflectance of the short infrared band of the Sentinel-2 satellite remote sensing image, and MIR is the surface reflectance of the mid-infrared band of the Sentinel-2 satellite remote sensing image.

[0020] Furthermore, the temporary land use range on the Sentinel-2 satellite remote sensing image includes multiple pixels, and the NDVI, NDBI, and NDSI values ​​calculated within the vector range are averaged;

[0021] For the time series NDVI, NDBI, and NDSI remote sensing characteristic curves generated for each plot, the maximum and minimum normalized vegetation index MMNDVI, the maximum and minimum normalized building index MMNDBI, and the maximum and minimum normalized soil index MMNDSI of the plot in each phase are obtained respectively by using the maximum and minimum normalization;

[0022] The expression of MMNDVI characteristic index is as follows:

[0023]

[0024] In the formula, MMNDVI represents the maximum and minimum normalized vegetation index, NDVI represents the normalized vegetation index at a certain moment; max1 represents the maximum value in the NDVI time series data, and min1 represents the minimum value in the NDVI time series data; through the maximum and minimum normalization processing, the values ​​of all calculated MMNDVI remote sensing characteristic indexes are between 0 and 1;

[0025] The expression of MMNDBI characteristic index is as follows:

[0026]

[0027] In the formula, MMNDBI represents the maximum and minimum normalized building index, NDBI represents the normalized building index at a certain moment; max2 represents the maximum value in the NDBI time series data, and min2 represents the minimum value in the NDBI time series data; through the maximum and minimum normalization processing, the values ​​of all calculated MMNDBI remote sensing characteristic indexes are between 0 and 1;

[0028] The MMNDSI characteristic index expression is as follows:

[0029]

[0030] In the formula, MMNDSI represents the maximum and minimum normalized soil index, NDSI represents the normalized soil index at a certain moment, max3 represents the maximum value in the NDSI time series data, and min3 represents the minimum value in the NDSI time series data. Through the maximum and minimum normalization processing, the values ​​of all calculated MMNDSI remote sensing characteristic indexes are between 0 and 1.

[0031] According to the determined temporary land vector, the usage cycle of temporary land and the calculated NDVI, NDBI, and NDSI remote sensing characteristic indexes of each patch, the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of the temporary land category were drawn respectively through mean processing and maximum and minimum normalization processing, and the typical time series sample characteristics of temporary land were obtained; 30 sample points of non-temporary land were randomly selected to draw the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of non-temporary land throughout the year.

[0032] In step S2, judging whether the land parcel is temporary land based on the DTW distance criterion includes:

[0033] Calculate the three remote sensing characteristic indexes including NDVI, NDBI and NDSI of the plot to be judged in each phase, and generate the time series characteristic curve of the three remote sensing characteristic indexes including MMNDVI, MMNDBI and MMNDSI of the land patch to be judged;

[0034] By comparing the DTW distance similarity between the time series characteristic curves of the three remote sensing characteristic indexes, including MMNDVI, MMNDBI and MMNDSI, of the plot to be judged and the time series characteristic curves of typical temporary land and non-temporary land, it is determined whether the plot to be judged is temporary land.

[0035] Furthermore, after obtaining the typical time series sample characteristics of non-temporary land from the time series characteristic curves of the three remote sensing characteristic indexes, the dynamic time warping (DTW) method is used to calculate the similarity of the time series characteristic curves of the three remote sensing characteristic indexes to determine whether the undetermined land parcel is temporary land.

[0036] For two time series X n ={c1,c2…c n} and Y m ={q1,q2…q m}, the lengths are |X| and Y| respectively, calculate the similarity of the two time series, and determine the regular path W=W1,W2…W K ; Among them, Max(|X|,|Y|)≤K≤|X|+|Y|, W K The form is (i, j), i represents the i coordinate in X, and j represents the j coordinate in Y; the normalization path W starts from W1 = (1,1) and ends at W1 = (|X|, |Y|) to ensure that every coordinate in X and Y appears in W; where c n For the first time series X n The value of the nth phase in m For the second time series Y m The value of the mth phase in K is the Kth point in the path, Max(|X|,|Y|) is the larger value of the two time series lengths, and |X|+|Y| is the sum of the two time series lengths;

[0037] In the regularized path W, i and j of W(i,j) are monotonically increasing, expressed as:

[0038] W K =(i,j)

[0039] W K+1 =(i′,j′)

[0040] i≤i′≤i+1,j≤j′≤j+1

[0041] Where i and j represent the position index of the Kth point in the first time series and the second time series respectively; W K+1 is the K+1th point in the path, i' and j' represent the position index of the K+1th point in the first time series and the second time series respectively, i≤i'≤i+1,j≤j'≤j+1 means that the movement of each step in the path is monotonically increasing;

[0042] The shortest optimal normalization path is the DTW distance, which is:

[0043] D(i,j)=Dist(i,j)+min[D(i-1,j),D(i,j-1),D(i-1,j-1)]

[0044] Where D(i,j) is the accumulated DTW distance between the time series at indexes i and j, Dist(i,j) is the local Euclidean distance between the i-th point in the first time series and the j-th point in the second time series, min[] represents the smallest of the three possible paths, i.e., the minimum accumulated distance from the top, left, or diagonal direction, D(i-1,j) is the accumulated DTW distance between the time series at indexes i-1 and j, D(i,j-1) is the accumulated DTW distance between the time series at indexes i and j-1, and D(i-1,j-1) is the accumulated DTW distance between the time series at indexes i-1 and j-1.

[0045] The DTW distances between the time series characteristic curves of the three remote sensing characteristic indexes of the pending land parcel and the time series characteristic curves of the three remote sensing characteristic indexes of the temporary land and the non-temporary land are calculated respectively. When the DTW distance with the temporary land is smaller, the pending land parcel is judged as temporary land, otherwise it is judged as non-temporary land.

[0046] In step S3, it is determined whether the temporary land has been reclaimed, including:

[0047] According to the time series characteristic curves of three remote sensing characteristic indices including MMNDVI, MMNDBI and MMNDSI of temporary land, the slopes of the three remote sensing characteristic indices at the expiration of the temporary land use period and the expiration of the reclamation period are calculated respectively; among which, the expiration of the reclamation period refers to the expiration of one year of the use period; the difference in the slopes of the three remote sensing characteristic indices is combined to determine whether the temporary land with the expiration of the reclamation period has completed land reclamation.

[0048] Furthermore, the relevant information for determining whether the temporary land with expiration of reclamation period has completed land reclamation includes: for the patches of temporary land with expiration of reclamation period, the three characteristic curve slopes are calculated using the three remote sensing characteristic indexes at the expiration time of the time series characteristic curve and the expiration time of reclamation period, respectively. The calculation expression is as follows:

[0049]

[0050] In the formula, K MMNDVI ,K MMNDBI ,K MMNDSI are the curvatures of the time series characteristic curves of the remote sensing characteristic indexes MMNDVI, MMNDBI, and MMNDSI; MMNDVI n ,MMNDBI n ,MMNDSI n are the values ​​of MMNDVI, MMNDBI, and MMNDSI remote sensing characteristic indexes on the time series characteristic curve at the time of reclamation expiration. m ,MMNDBI m,MMNDSI m are the values ​​of the time of expiration of the use period on the time series characteristic curves of the remote sensing characteristic indexes of MMNDVI, MMNDBI, and MMNDSI, respectively; m and m are the time of expiration of the reclamation period and the time of expiration of the use period, respectively;

[0051] The slope values ​​of the time series characteristic curves of three remote sensing characteristic indices, including MMNDVI, MMNDBI and MMNDSI, are combined to determine whether land reclamation is completed;

[0052] When K is satisfied MMNDVI >δ1,K MMNDBI <δ2,K MMNDSI <δ3 indicates that the land has been reclaimed and the reclamation effect is good;

[0053] When K is satisfied MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI ≥δ3 indicates that the land reclamation has not been completed on the plot;

[0054] When K cannot be satisfied at the same time MMNDVI >δ1,K MMNDBI <δ2,K MMNDSI <δ3 or K MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI When ≥δ3, manual interpretation is performed;

[0055] Among them, δ1, δ2, and δ3 are the thresholds of the slopes of the three remote sensing characteristic indexes used to determine whether land reclamation has been completed, which are calculated based on 30 temporary land samples that have completed land reclamation. In practical applications, they are obtained by statistical calculation based on temporary land samples in different regions.

[0056] Another object of the present invention is to provide a temporary land use interpretation and reclamation verification system, which is used to control the temporary land use interpretation and reclamation verification method, and the system includes:

[0057] The time series remote sensing index sample feature generation module is used to calculate and generate sample features of typical time series characteristic curves of remote sensing characteristic index of temporary land and non-temporary land by using GEE cloud computing platform and Sentinel-2 satellite remote sensing images;

[0058] The temporary land use judgment module is used to calculate the time series characteristic curve of the remote sensing characteristic index of the land to be judged, and after linear interpolation or linear compression, calculate the DTW similarity distance between the time series characteristic curve of the remote sensing characteristic index of the land to be judged and the typical time series characteristic curves of temporary land and non-temporary land, and judge whether the land is temporary land based on the DTW distance criterion;

[0059] The land reclamation judgment module is used to calculate the slope of the time series characteristic curve for temporary land that has expired, based on the remote sensing characteristic index at the time of expiration of the use period and the time of expiration of the reclamation period, to determine whether the temporary land has completed land reclamation.

[0060] Existing land reclamation verification schemes for temporary land are mainly based on on-site verification, and the verification methods directly based on satellite remote sensing images or drone remote sensing are also mainly judged by manual visual interpretation, which requires a lot of manpower and time costs. In response to the above problems, the advantages and positive effects of the present invention are as follows: the present invention analyzes the difference in the time series characteristic curves of the remote sensing characteristic index of temporary land and non-temporary land, and uses the difference in the characteristic curves of the remote sensing characteristic index of temporary land that has completed land reclamation and temporary land that has not been reclaimed on the time series satellite remote sensing images. Rapid screening of temporary land and rapid, automated judgment of whether land reclamation has been completed have been achieved, breaking through the current technical bottleneck of relying on manual visual interpretation when verifying temporary land directly based on satellite remote sensing images. At the same time, this method also effectively avoids the problem of subjective verification results caused by human subjective inference.

[0061] The satellite remote sensing images used in this method are Sentinel-2 optical remote sensing images, and can also be expanded to other optical satellite remote sensing data with red bands, near infrared bands, short infrared bands, and mid-infrared bands. In addition, due to the limitations of the spatial resolution of Sentinel-2 optical satellite remote sensing images and the scale of temporary land use, this method is more suitable for the identification of temporary land use with an area greater than 0.7 hectares and less than 8 hectares and automatic verification of land reclamation.

[0062] The present invention was verified by an example of 50 temporary lands that have been confirmed to have completed land reclamation in Zhejiang Province. The accuracy rate of temporary land judgment reached 84%, and the error rate was only 4%. The accuracy rate of temporary land judgment for land reclamation reached 86%. It is proved that the method of the present invention can effectively identify temporary land and judge whether the temporary land has completed land reclamation within the prescribed period, which greatly reduces the manpower cost. The present invention is a method for indirectly using satellite remote sensing images to screen whether temporary land has completed land reclamation, which provides a new technical method for monitoring and supervising temporary land.

[0063] For temporary land that occupies cultivated land, forest land, grassland, agricultural land, etc., the present invention provides a new method flow to simplify the method of verifying whether the temporary land patches have completed land reclamation, which can greatly reduce the time cost of manual visual interpretation and verification of land reclamation after the expiration of the temporary land use period. The temporary land interpretation and reclamation verification method based on time-series satellite remote sensing images of the present invention uses time-series satellite remote sensing images to continuously track and monitor the problems existing in the use of temporary land, especially to grasp the land reclamation situation after the expiration of the temporary land use period, which can effectively supplement the coverage and timeliness of the traditional supervision model, save a lot of manpower and material resources, and significantly improve the government's supervision level and decision-making analysis capabilities.

[0064] The present invention realizes the automation and process of temporary land verification by using remote sensing images and related algorithms, which can significantly improve the efficiency and accuracy of temporary land verification, and can obtain the results of most temporary land interpretation and reclamation verification without manual verification, significantly reducing labor costs. It is expected to bring significant economic benefits and social value, especially in the fields of land resource management, environmental protection, etc., and has broad commercial application prospects.

[0065] The present invention realizes for the first time the automation of temporary land interpretation and land reclamation verification based on time-series remote sensing images, and integrates a temporary land verification system, filling the gap in the automation of temporary land supervision and verification based on time-series remote sensing images.

[0066] Traditional temporary land verification relies on manual on-site investigation, which is time-consuming, labor-intensive and easily affected by human factors. By using time-series remote sensing images and related algorithms, the present invention builds a temporary land interpretation and reclamation verification system based on the GEE cloud platform, realizing the automation of temporary land interpretation and land reclamation verification process, significantly reducing labor costs and time costs, and overcoming the limitations of traditional temporary land verification technical means. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] The accompanying drawings herein are incorporated in and constitute a part of the specification, illustrate embodiments consistent with the present disclosure, and together with the description, serve to explain the principles of the present disclosure;

[0068] Figure 1 It is a flow chart of a temporary land use interpretation and reclamation verification method provided by an embodiment of the present invention;

[0069] Figure 2 It is a time series characteristic curve diagram of MMNDVI, MMNDBI, and MMNDSI remote sensing characteristic indexes of temporary land use throughout the entire period provided by an embodiment of the present invention;

[0070] Figure 3It is a full-year time series characteristic curve diagram of MMNDVI, MMNDBI, and MMNDSI remote sensing characteristic indexes of non-temporary land (vegetation area) provided by an embodiment of the present invention;

[0071] Figure 4 It is a flow chart of a method for distinguishing temporary land use according to a time series characteristic curve combining three remote sensing characteristic indexes provided by an embodiment of the present invention;

[0072] Figure 5 It is a flow chart of a method provided by an embodiment of the present invention for determining whether land reclamation of temporary land with an expired reclamation period has been completed based on the slope of a time series characteristic curve combining three remote sensing characteristic indexes. DETAILED DESCRIPTION

[0073] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention, so the present invention is not limited by the specific implementation disclosed below.

[0074] The innovation of the present invention lies in that: the present invention utilizes time-series Sentinel-2 optical remote sensing images and related algorithms, and for the first time proposes to combine MMNDVI, MMNDBI, MMNDSI and other three remote sensing characteristic indexes and the slope value of the time-series characteristic curve to automatically judge whether the temporary land has completed land reclamation. Based on the GEE cloud computing platform, the present invention realizes the automation of temporary land interpretation and land reclamation verification for the first time, overcoming the problem that the remote sensing verification of temporary land used to rely on manual interpretation in the past.

[0075] The present invention is the first to jointly apply remote sensing characteristic indices calculated using three formulas, including NDVI, NDBI and NDSI, to the feature construction and judgment of temporary use, and uses these as criteria for judging temporary land use, thus solving the problem that temporary land use can only rely on manual interpretation.

[0076] In addition, since the temporary land use range on the Sentinel-2 satellite remote sensing image generally includes multiple pixels, it is necessary to average the NDVI, NDBI, and NDSI values ​​calculated within the vector range. In order to weaken the impact of different geographical locations on the remote sensing index, for the time series NDVI, NDBI, and NDSI remote sensing characteristic curves generated for each plot, the maximum and minimum normalized vegetation index MMNDVI, the maximum and minimum normalized building index MMNDBI, and the maximum and minimum normalized soil index MMNDSI of the plot in each phase are obtained using maximum and minimum normalization.

[0077] Embodiment 1, as Figure 1 As shown, the temporary land use interpretation and reclamation verification method provided by the embodiment of the present invention includes the following steps:

[0078] S1. Using the GEE cloud computing platform and Sentinel-2 satellite remote sensing images, calculate and generate sample characteristics of typical time series characteristic curves of remote sensing characteristic indexes of temporary land and non-temporary land;

[0079] According to the vector range shp file or kml file of temporary land, the GEE (Google Earth Engine, GEE) cloud platform and Sentinel-2 satellite remote sensing images were used to obtain the time series NDVI, NDBI and NDSI remote sensing characteristic indices of temporary land within the land use cycle and the full year time series NDVI, NDBI and NDSI remote sensing characteristic indices of non-temporary land. The time series MMNDVI, MMNDBI and MMNDSI values ​​of typical temporary and non-temporary land patches were obtained by performing maximum and minimum normalization processing on the time series characteristic curves. The sample characteristics of the typical time series characteristic curves of the three remote sensing characteristic indices, MMNDVI, MMNDBI and MMNDSI, of temporary and non-temporary land categories were then obtained by mean calculation.

[0080] S2, calculating the time series characteristic curve of the remote sensing characteristic index of the plot to be judged, and after linear interpolation or linear compression, calculating the DTW similarity distance between the time series characteristic curve of the remote sensing characteristic index of the plot to be judged and the typical time series characteristic curves of temporary land and non-temporary land, and judging whether the plot is temporary land based on the DTW distance criterion;

[0081] To determine whether the pending land parcel is temporary land, three remote sensing characteristic indexes, namely NDVI, NDBI and NDSI, are calculated for each time phase of the pending land parcel, and time series characteristic curves of three remote sensing characteristic indexes, namely MMNDVI, MMNDBI and MMNDSI, are generated for the land patch to be determined. By comparing the DTW distance similarity between the time series characteristic curves of the three remote sensing characteristic indexes of the pending land parcel and the time series characteristic curves of typical temporary land and non-temporary land, it is determined whether the pending land parcel is temporary land.

[0082] S3. For temporary land that has expired, the slope of the time series characteristic curve is calculated based on the remote sensing characteristic index at the time of expiration of the use period and the time of expiration of the reclamation period to determine whether the temporary land has completed land reclamation;

[0083] Determine whether the temporary land with expiration of reclamation period has completed land reclamation. According to the time series characteristic curves of the three remote sensing characteristic indexes of temporary land, MMNDVI, MMNDBI and MMNDSI, the slopes of the three remote sensing characteristic indexes at the time of expiration of the temporary land use period and the time of expiration of reclamation period (generally refers to the expiration of one year of use) are calculated respectively. Combined with the difference of the slopes of the three remote sensing characteristic indexes, it is possible to determine whether the temporary land with expiration of reclamation period has completed land reclamation and other related information.

[0084] Before step S1, it is also necessary to perform preprocessing work such as image screening, cloud removal, and image cropping on the Sentinel-2 satellite remote sensing images to obtain temporary and non-temporary land samples with high availability and generate time series information on the entire life cycle of temporary land. The screening of Sentinel-2 satellite remote sensing images is based on the GEE platform. The date is set according to the temporary land vector range, approval time, and expiration time. Image data with cloud cover less than 20% during the temporary land use period and land reclamation period are screened for subsequent calculations. Then, the obtained Sentinel-2 satellite remote sensing images are declouded using the QA60 band to reduce the impact of cloud cover on data statistics.

[0085] In step S1, remote sensing characteristic indexes such as NDVI, NDBI and NDSI of the Sentinel-2 satellite remote sensing image of each phase are calculated. The Sentinel-2 satellite sensor includes 12 spectral bands. The normalized vegetation index NDVI (Normalized DifferentialVegetation Index, NDVI), the normalized building index NDBI (Normalized Differential Built Up Index, NDBI) and the normalized soil index NDSI (Normalized Differential Soil Index, NDSI) are calculated respectively by band operation to obtain the NDVI, NDBI and NDSI values ​​of satellite remote sensing data in different phases.

[0086] Among them, the calculation expressions of NDVI, NDBI and NDSI are as follows:

[0087]

[0088] Where NIR is the surface reflectance of the near infrared band of the Sentinel-2 satellite remote sensing image, RED is the surface reflectance of the red band of the Sentinel-2 satellite remote sensing image, SWIR is the surface reflectance of the short infrared band of the Sentinel-2 satellite remote sensing image, and MIR is the surface reflectance of the mid-infrared band of the Sentinel-2 satellite remote sensing image.

[0089] In addition, since the temporary land use range on the Sentinel-2 satellite remote sensing image generally includes multiple pixels, it is necessary to average the NDVI, NDBI, and NDSI values ​​calculated within the vector range.

[0090] For the time series NDVI, NDBI, and NDSI remote sensing characteristic curves generated for each plot, the maximum and minimum normalized vegetation index MMNDVI, the maximum and minimum normalized building index MMNDBI, and the maximum and minimum normalized soil index MMNDSI of the plot in each phase were obtained respectively using maximum and minimum normalization.

[0091] MMNDVI (Max-Min Normalization Differential Vegetation Index, MMNDVI), called the maximum and minimum normalized vegetation index, is a time-series modification of the NDVI index. It is a remote sensing characteristic index that uses index volatility to measure changes in surface vegetation coverage. It can objectively reflect changes in vegetation coverage, weaken the impact of different geographical locations on NDVI, and only rely on the law of index changes to judge changes in surface vegetation within a certain period of time. Generally speaking, the higher the MMNDVI value, the richer the surface vegetation. The corresponding MMNDVI characteristic index expression is as follows:

[0092]

[0093] In the formula, MMNDVI represents the maximum and minimum normalized vegetation index, NDVI represents the normalized vegetation index at a certain moment; max1 represents the maximum value in the NDVI time series data, and min1 represents the minimum value in the NDVI time series data; through the maximum and minimum normalization processing, the values ​​of all calculated MMNDVI remote sensing characteristic indexes are between 0 and 1;

[0094] MMNDBI (Max-Min Normalization Differential Built Up Index, MMNDBI), called the maximum and minimum normalized building index, is a time-series modification of the NDBI index. It is a remote sensing characteristic index used to measure the density of surface buildings. It reflects the changes in the density or richness of surface buildings, weakens the impact of different geographical locations on NDBI, and only relies on the index change law to judge the change information of surface buildings in a certain period of time. Generally speaking, the higher the MMNDBI value, the richer the surface buildings. The corresponding MMNDBI characteristic index expression is as follows:

[0095]

[0096] In the formula, MMNDBI represents the maximum and minimum normalized building index, NDBI represents the normalized building index at a certain moment; max2 represents the maximum value in the NDBI time series data, and min2 represents the minimum value in the NDBI time series data; through the maximum and minimum normalization processing, the values ​​of all calculated MMNDBI remote sensing characteristic indexes are between 0 and 1;

[0097] MMNDSI (Max-Min Normalization Differential Soil Index, MMNDSI), called the maximum and minimum normalized soil index, is a remote sensing characteristic index used to measure the abundance of bare soil on the surface. It is a time-series modification of the NDSI index, which weakens the impact of different geographical locations on the NDSI and only relies on the index change law to judge the change information of bare soil in a certain period of time. Generally speaking, the higher the MMNDSI value, the higher the bare soil content on the surface. The corresponding MMNDSI characteristic index expression is as follows:

[0098]

[0099] In the formula, MMNDSI represents the maximum and minimum normalized soil index, NDSI represents the normalized soil index at a certain moment, max3 represents the maximum value in the NDSI time series data, and min3 represents the minimum value in the NDSI time series data. Through the maximum and minimum normalization processing, the values ​​of all calculated MMNDSI remote sensing characteristic indexes are between 0 and 1.

[0100] According to the determined temporary land vector, the use cycle of temporary land, and the calculated NDVI, NDBI, and NDSI remote sensing characteristic indexes of each patch, the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of the temporary land category are drawn respectively through mean processing and maximum and minimum normalization processing, so as to obtain the typical time series sample characteristics of temporary land; at the same time, 30 sample points of non-temporary land are randomly selected, and the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of non-temporary land throughout the year are drawn according to the above method.

[0101] In step S2, according to the above method, the time series characteristic curves of the three remote sensing characteristic indexes of the undetermined land parcel are generated, that is, the typical time series sample characteristics of non-temporary land are obtained; then, the time series characteristic curves of the three remote sensing characteristic indexes are combined, and the dynamic time warping DTW (Dynamic TimeWarping, DTW) method is used to calculate the similarity to determine whether the undetermined land parcel is temporary land. DTW distance judgment similarity can be used for data of unequal length and focuses on feature matching. Compared with Euclidean distance, it is more natural to align data.

[0102] DTW uses a time warping function W(n) that satisfies certain conditions to describe the time correspondence between the test template and the reference template, and solves the warping function corresponding to the minimum cumulative distance when the test template and the reference template are matched. n ={c1,c2…c n} and Y m ={q1,q2…q m}, the lengths are X| and |Y| respectively, calculate the similarity of the two time series, and determine the regular path W=W1,W2…W K ; Among them, Max(|X|,|X|)≤K≤|X|+|Y|, W K The form is (i, j), i represents the i coordinate in X, and j represents the j coordinate in Y; the normalization path W starts from W1 = (1,1) and ends at W1 = (|X|, |Y|) to ensure that every coordinate in X and Y appears in W; where c n For the first time series X n The value of the nth phase in m For the second time series Y m The value of the mth phase in K is the Kth point in the path, Max(|X|,|Y|) is the larger value of the two time series lengths, and |Y|+|Y| is the sum of the two time series lengths;

[0103] In the regularized path W, i and j of W(i,j) are monotonically increasing, expressed as:

[0104] W K =(i,j)

[0105] W K+1 =(i′,j′)

[0106] i≤i′≤i+1,j≤j′≤j+1

[0107] Where i and j represent the position index of the Kth point in the first time series and the second time series respectively; W K+1 is the K+1th point in the path, i' and j' represent the position index of the K+1th point in the first time series and the second time series respectively, i≤i'≤i+1,j≤j'≤j+1 means that the movement of each step in the path is monotonically increasing;

[0108] The shortest optimal normalization path is the DTW distance, which is:

[0109] D(i,j)=Dist(i,j)+min[D(i-1,j),D(i,j-1),D(i-1,j-1)]

[0110] Where D(i,j) is the accumulated DTW distance between the time series at index i and j, Dist(i,j) is the local Euclidean distance between the i-th point in the first time series and the j-th point in the second time series, min[] represents the smallest one among the three possible paths, that is, the minimum accumulated distance from the top, left or diagonal direction, D(i-1,j) is the accumulated DTW distance between the time series at index i-1 and j, D(i,j-1) is the accumulated DTW distance between the time series at index i and j-1, and D(i-1,j-1) is the accumulated DTW distance between the time series at index i-1 and j-1;

[0111] The DTW distances between the time series characteristic curves of the three remote sensing characteristic indexes of the pending land parcel and the time series characteristic curves of the three remote sensing characteristic indexes of the temporary land and the non-temporary land are calculated respectively. When the DTW distance with the temporary land is smaller, the pending land parcel is judged as temporary land, otherwise it is judged as non-temporary land.

[0112] In step S3, for the patch of temporary land determined to have expired reclamation period, the three characteristic curve slopes are calculated using the three remote sensing characteristic indices at the time of expiration of the use period and the time of expiration of the reclamation period of the time series characteristic curve, respectively. The calculation expression is as follows:

[0113]

[0114] In the formula, K MMNDVI ,K MMNDBI ,K MMNDSI are the curvatures of the time series characteristic curves of the remote sensing characteristic indexes MMNDVI, MMNDBI, and MMNDSI; MMNDVI n ,MMNDBI n ,MMNDSI n are the values ​​of MMNDVI, MMNDBI, and MMNDSI remote sensing characteristic indexes on the time series characteristic curve at the time of reclamation expiration. m ,MMNDBI m ,MMNDSI m are the values ​​of the time of expiration of the use period on the time series characteristic curves of the remote sensing characteristic indexes of MMNDVI, MMNDBI and MMNDSI, respectively; n and m are the time of expiration of the reclamation period and the time of expiration of the use period, respectively;

[0115] The slope values ​​of the time series characteristic curves of three remote sensing characteristic indices, including MMNDVI, MMNDBI and MMNDSI, are combined to determine whether land reclamation is completed;

[0116] When K is satisfied MMNDVI >δ1,K MMNDBI <δ2,KMMNDSI <δ3 indicates that the land has been reclaimed and the reclamation effect is good;

[0117] When K is satisfied MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI ≥δ3 indicates that the land reclamation has not been completed on the plot;

[0118] When K cannot be satisfied at the same time MMNDVI >δ1,K MMNDBI <δ2,K MMNDSI <δ3 or K MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI When ≥δ3, manual interpretation is performed;

[0119] Among them, δ1, δ2, and δ3 are the thresholds of the slopes of the three remote sensing characteristic indexes used to determine whether land reclamation has been completed, which are calculated based on 30 temporary land samples that have completed land reclamation. In practical applications, they are obtained by statistical calculation based on temporary land samples in different regions.

[0120] Embodiment 2, the method for identifying temporary land use and verifying land reclamation based on time-series satellite remote sensing images provided by the embodiment of the present invention specifically includes the following steps 1 to 4:

[0121] Step 1, the present invention first calculates the MMNDVI, MMNDBI and MMNDSI remote sensing characteristic indexes of the temporary land category, which are used to extract the fluctuation characteristics of the time series characteristic curve and weaken the interference of regional position differences, so as to obtain the typical time series characteristics of the temporary land, which specifically includes the following steps:

[0122] (1.1) Based on the Sentinel-2 optical remote sensing image patches of temporary land use in Zhejiang Province, data integration is performed on the patches according to approval time and expiration time.

[0123] (1.2) Manually check temporary map areas to ensure that the map information is accurate.

[0124] (1.3) The remote sensing image patches of 10 temporary lands that have been confirmed to have completed land reclamation through field investigation were selected for feature statistics.

[0125] (1.4) Load temporary land use patches one by one on the Google Earth Engine (GEE) platform, and use the ee.ImageCollection function to filter images with cloud cover less than 20% based on the approval and expiration time of temporary land use.

[0126] (1.5) The QA60 band of the GEE platform is used to remove the cloud from the Sentinel-2 satellite remote sensing images of each phase.

[0127] (1.6) Calculate the MMNDVI, MMNDBI and MMNDSI remote sensing characteristic indexes of temporary land patches, and draw the time series characteristic curves of NDVI, NDBI, NDSI, etc. from the approval of temporary land use to the reclamation acceptance stage.

[0128] For Sentinel-2 optical remote sensing satellite data, the MMNDVI index is calculated using the 8th and 4th bands of Sentinel-2 optical remote sensing satellite data and the NDVI time series characteristic curve. The calculation formula is:

[0129]

[0130] Where B8 and B4 are the near-infrared band and red light band of Sentinel-2 optical remote sensing satellite data, max1 represents the maximum value in the NDVI time series data, and min1 represents the minimum value in the NDVI time series data.

[0131] MMNDBI uses the 11th and 8th bands of Sentinel-2 optical remote sensing satellite data and the NDBI time series characteristic curve for calculation. The calculation formula is:

[0132]

[0133] Where B11 and B18 are the short infrared band and near infrared band of Sentinel-2 optical remote sensing satellite data, max2 represents the maximum value in the NDBI time series data, and min2 represents the minimum value in the NDBI time series data.

[0134] MMNDSI is calculated using the 12th and 8th bands of the Sentinel-2 optical remote sensing satellite data and the NDSI time series characteristic curve. The calculation formula is:

[0135]

[0136] Where B12 and B8 are the mid-infrared band and near-infrared band of Sentinel-2 optical remote sensing satellite data, max3 represents the maximum value in the NDSI time series data, and min3 represents the minimum value in the NDSI time series data.

[0137] (1.7) Perform linear interpolation and linear compression on the time series characteristic curves of the 10 temporary land patches that have been determined, so that the average length of their horizontal axes is 50.

[0138] (1.8) Taking the characteristic mean of 10 temporary land spots at each horizontal axis position, we can obtain the time series characteristic curves of the three remote sensing characteristic indices of temporary land use, namely, MMNDVI, MMNDBI, and MMNDSI, as shown in the figure: Figure 2 shown.

[0139] Next, the present invention uses the GEE cloud platform to calculate the time series characteristic curves of three remote sensing characteristic indexes, namely MMNDVI, MMNDBI and MMNDSI, for non-temporary land (vegetation) throughout the year.

[0140] (1.9) In the GEE cloud platform, 30 non-temporary land (vegetation) samples were randomly selected in the vicinity of the identified temporary land.

[0141] (1.10) Taking the 2019 Sentinel-2 remote sensing data as an example, the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of 30 non-temporary land use samples are calculated respectively.

[0142] (1.11) The time series characteristic curves of the 30 non-temporary land use samples are linearly compressed so that the length of the horizontal axis is 12 (i.e., one value per month).

[0143] (1.12) Take the characteristic mean of 30 non-temporary land samples at each horizontal axis position to obtain the annual time series characteristic curves of MMNDVI, MMNDBI and MMNDSI of non-temporary land categories, as shown in Figure 3 shown.

[0144] Step 2, such as Figure 4 As shown in the figure, for each patch to be judged, the above method is used to generate the time series characteristic curve of the three remote sensing characteristic indexes of the plot, and then the DTW distance algorithm is used to calculate the similarity between the plot to be judged and the temporary land and non-temporary land, so as to judge whether the patch is a temporary land based on the similarity. The specific steps are as follows:

[0145] (2.1) First, the time range of the MMNDVI, MMNDBI and MMNDSI time series characteristic curves of the land to be judged should be determined based on the approval start time and the reclamation acceptance time.

[0146] (2.2) Input the plot to be judged into the GEE cloud platform and obtain the MMNDVI, MMNDBI and MMNDSI time series characteristic curves of the plot to be judged.

[0147] (2.3) Calculate the similarity between the MMNDVI time series characteristic curve of the plot to be judged, the MMNDVI time series characteristic curve of the temporary land, and the MMNDVI time series characteristic curve of the non-temporary land, that is, the DTW distance, and obtain D 临时用地MMNDVI and D 非临时用地MMNDVICalculate the similarity between the MMNDBI time series characteristic curve of the segment to be judged, the MMNDBI time series characteristic curve of the temporary land, and the MMNDBI time series characteristic curve of the non-temporary land, that is, the DTW distance, and obtain D 临时用地MMNDBI and D 非临时用地MMNDBI Calculate the similarity between the MMNDSI timing characteristic curve of the segment to be judged and the MMNDSI timing characteristic curve of the temporary land and the MMNDSI timing characteristic curve of the non-temporary land, that is, the DTW distance, and obtain D 临时用地MMNDSI and D 非临时用地MMNDSI .

[0148] (2.4) Based on the calculated 6 DTW distance results, the plot to be judged is automatically judged. The judgment rules are as follows: Figure 4 As shown: For condition D 临时用地MMNDVI <D 非临时用地MMNDVI , D 临时用地MMNDBI <D 非临时用地MMNDBI , D 临时用地MMNDSI <D 非临时用地MMNDSI , when two of them are met, it is judged as temporary land use; when D 临时用地MMNDVI >D 非临时用地MMNDVI , D 临时用地MMNDBI >D 非临时用地MMNDBI , D 临时用地MMNDSI >D 非临时用地MMNDSI时 , it is judged as non-temporary land use; other situations need to be transferred to manual interpretation.

[0149] (2.5) Taking the temporary land patches in Zhejiang Province as the data source, 50 plots were selected as verification data, and the accuracy and error rate of the discrimination method were calculated. The accuracy rate of the automatic discrimination of temporary land in the present invention reached 84%, and the error rate was only 4%.

[0150] Step 3, such as Figure 5 As shown in the figure, the characteristic morphology of the MMNDVI, MMNDBI and MMNDSI time series characteristic curves of the temporary land patches is used to determine whether the land reclamation is completed. The specific steps are as follows:

[0151] (3.1) For temporary land that has expired, extract the MMNDVI, MMNDBI and MMNDSI characteristic values ​​of the temporary land patch at the time of expiration and the MMNDVI, MMNDBI and MMNDSI characteristic values ​​at the time of expiration of the reclamation period.

[0152] (3.2) Calculate the slopes of MMNDVI, MMNDBI and MMNDSI, which are K MMNDVI ,K MMNDBI ,K MMNDSI , the calculation expression is as follows:

[0153]

[0154] In the formula, K MMNDVI ,K MMNDBI ,K MMNDSI are the curvatures of the time series characteristic curves of the remote sensing characteristic indexes MMNDVI, MMNDBI, and MMNDSI; MMNDVI n ,MMNDBI n ,MMNDSI n are the values ​​of MMNDVI, MMNDBI, and MMNDSI remote sensing characteristic indexes on the time series characteristic curve at the time of reclamation expiration. m ,MMNDBI m ,MMNDSI m are the values ​​of the time of expiration of the use period on the time series characteristic curves of the remote sensing characteristic indexes of MMNDVI, MMNDBI and MMNDSI, respectively; n and m are the time of expiration of the reclamation period and the time of expiration of the use period, respectively;

[0155] (3.3) The three slope values ​​are used to determine whether land reclamation is complete. The judgment rules are as follows: Figure 5 As shown: When K MMNDVI >δ1,K MMNDBI <δ2,K MMNDSI <δ3, indicating that the land reclamation of the plot is well completed; when K MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI ≥δ3, indicating that the land reclamation of the plot has not been completed; when only one of the conditions is met, manual interpretation is required.

[0156] (3.4) Taking the temporary land patches in Zhejiang Province as the data source, 50 temporary land patches that have completed land reclamation are selected as verification data, and the accuracy of automatic judgment of land reclamation is calculated. The accuracy of the automatic judgment of whether land reclamation is completed by the present invention reaches 86%.

[0157] Step 4: Result evaluation.

[0158] Currently, the existing temporary land reclamation verification adopts the method of on-site verification or manual visual interpretation of high-resolution satellite images, and there is no automated temporary land identification and temporary land reclamation identification technology. The present invention proposes a new temporary land identification and land reclamation technology process.

[0159] (4.1) Qualitative evaluation: Currently, manual visual interpretation is used, and there is no automated technology for temporary land identification and temporary land reclamation identification. Based on the time-series Sentinel-2 optical satellite remote sensing images, this invention calculates three remote sensing characteristic indexes, namely MMNDVI, MMNDBI and MMNDSI, constructs a time-series characteristic curve, and combines the DTW algorithm to achieve automatic verification of temporary land identification and land reclamation, overcoming the drawbacks of relying solely on manual verification in the past, which consumes a lot of manpower and time costs, and provides a new technical means for automatic monitoring and verification of land reclamation of temporary land.

[0160] (4.2) Quantitative evaluation: The accuracy of the automatic identification of temporary land of the present invention is 84%, which can meet the demand for automatic identification of temporary land. Although the accuracy of temporary land identification still needs to be improved, the error rate is only 4%. For temporary land that cannot be automatically identified by the present invention, manual verification can be continued, and the probability of direct error is low. The accuracy of the automatic verification of land reclamation of temporary land after the reclamation period expires of the present invention is 86%, with high accuracy and no verification errors. For temporary land that cannot automatically determine whether the land reclamation has been completed, manual verification can be continued.

[0161] Embodiment 3, the temporary land use interpretation and reclamation verification system provided by the embodiment of the present invention includes:

[0162] The time series remote sensing index sample feature generation module is used to calculate and generate sample features of typical time series characteristic curves of remote sensing characteristic index of temporary land and non-temporary land by using GEE cloud computing platform and Sentinel-2 satellite remote sensing images;

[0163] The temporary land use judgment module is used to calculate the time series characteristic curve of the remote sensing characteristic index of the land to be judged, and after linear interpolation or linear compression, calculate the DTW similarity distance between the time series characteristic curve of the remote sensing characteristic index of the land to be judged and the typical time series characteristic curves of temporary land and non-temporary land, and judge whether the land is temporary land based on the DTW distance criterion;

[0164] The land reclamation judgment module is used to calculate the slope of the time series characteristic curve for temporary land that has expired, based on the remote sensing characteristic index at the time of expiration of the use period and the time of expiration of the reclamation period, to determine whether the temporary land has completed land reclamation.

[0165] In order to further demonstrate the positive effects of the above embodiment, the present invention conducts the following experiments based on the above technical solution.

[0166] A temporary land patch in Zhejiang Province is taken as an example. The start and end time of the patch is from December 15, 2019 to December 29, 2021. The time series NDVI, NDBI and NDSI of the patch are obtained, and the time series feature length is 52; the time series features are normalized to the maximum and minimum to obtain MMNDVI, MMNDBI and MMNDSI, respectively, and linear compression is performed; then the DTW distance is calculated between the MMNDVI, MMNDBI and MMNDSI time series characteristic curves of the temporary land and the MMNDVI, MMNDBI and MMNDSI time series characteristic curves of the non-temporary land. After calculation, the DTW distances between the MMNDVI, MMNDBI and MMNDSI of the plot and those of the temporary land are 10.261, 9.828 and 13.757, respectively, and the DTW distances between the MMNDVI, MMNDBI and MMNDSI of the plot and those of the non-temporary land are 17.742, 25.102 and 18.183, respectively. 临时用地MMNDVI <D 非临时用地MMNDVI , D 临时用地MMNDBI <D 非临时用地MMNDBI , D 临时用地MMNDSI <D 非临时用地MMNDSI , the plot was interpreted as temporary land, which was correct. The slope of the characteristic curve was calculated using the MMNDVI, MMNDBI and MMNDSI values ​​at the beginning and end of the plot. Since the months were the same, δ1 = δ2 = δ3 = 0. After calculating K MMNDVI =-0.008≤δ1, K MMNDBI =0.001≥δ2, K MMNDSI =0.003≥δ3. It was determined that the land reclamation of the plot had not been completed. Combined with the actual situation, this determination was correct.

[0167] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by any technician familiar with the technical field within the technical scope disclosed by the present invention and within the spirit and principles of the present invention should be covered within the protection scope of the present invention.

Claims

1. A temporary land use interpretation and reclamation verification method, characterized in that: The method includes: S1. Using the GEE cloud computing platform and Sentinel-2 satellite remote sensing images, calculate and generate sample characteristics of typical time series characteristic curves of remote sensing characteristic indexes of temporary land and non-temporary land; S2, calculating the time series characteristic curve of the remote sensing characteristic index of the plot to be judged, and after linear interpolation or linear compression, calculating the DTW similarity distance between the time series characteristic curve of the remote sensing characteristic index of the plot to be judged and the typical time series characteristic curves of temporary land and non-temporary land, and judging whether the plot is temporary land based on the DTW distance criterion; S3. For temporary land that has expired, the slope of the time series characteristic curve is calculated based on the remote sensing characteristic index at the time of expiration of the use period and the time of expiration of the reclamation period to determine whether the temporary land has completed land reclamation.

2. The temporary land use interpretation and reclamation verification method according to claim 1 is characterized in that: Before step S1, the Sentinel-2 satellite remote sensing images are pre-processed by image screening, cloud removal and image cropping; temporary land and non-temporary land samples are obtained to generate time series information of the entire life cycle of temporary land; Sentinel-2 satellite remote sensing image screening is based on the GEE platform to set the date according to the temporary land vector range and the approval time and expiration time, and the image data with cloud cover less than 20% during the temporary land use period and the land reclamation period are screened, and the obtained Sentinel-2 satellite remote sensing images are de-clouded using the QA60 band.

3. The temporary land use interpretation and reclamation verification method according to claim 1 is characterized in that: In step S1, sample features of typical time series characteristic curves of remote sensing characteristic indexes of temporary land and non-temporary land are calculated and generated, including: According to the vector range shp file and kml file of temporary land, the GEE cloud platform and Sentinel-2 satellite remote sensing images are used to obtain the time series NDVI, NDBI and NDSI remote sensing characteristic indexes of temporary land during the land use period and the full year time series NDVI, NDBI and NDSI remote sensing characteristic indexes of non-temporary land. By performing maximum and minimum normalization processing on the time series characteristic curve, the time series MMNDVI, MMNDBI and MMNDSI values ​​of typical temporary land and non-temporary land patches were obtained; The sample characteristics of typical time series characteristic curves of three remote sensing characteristic indices, including MMNDVI, MMNDBI and MMNDSI, for temporary land use and non-temporary land use categories were obtained by mean calculation.

4. The temporary land use interpretation and reclamation verification method according to claim 3 is characterized in that: Obtain the time series NDVI, NDBI and NDSI remote sensing characteristic indexes of temporary land use during the land use period, including: The Sentinel-2 satellite sensor includes 12 spectral bands, and the normalized vegetation index NDVI, normalized building index NDBI and normalized soil index NDSI are calculated by band operation to obtain the NDVI, NDBI and NDSI values ​​of satellite remote sensing data at different phases. Among them, the calculation expressions of NDVI, NDBI and NDSI are as follows: Where NIR is the surface reflectance of the near infrared band of the Sentinel-2 satellite remote sensing image, RED is the surface reflectance of the red band of the Sentinel-2 satellite remote sensing image, SWIR is the surface reflectance of the short infrared band of the Sentinel-2 satellite remote sensing image, and MIR is the surface reflectance of the mid-infrared band of the Sentinel-2 satellite remote sensing image.

5. The temporary land use interpretation and reclamation verification method according to claim 4 is characterized in that: The temporary land use range on the Sentinel-2 satellite remote sensing image includes multiple pixel points, and the NDVI, NDBI, and NDSI values ​​calculated within the vector range are averaged; For the time series NDVI, NDBI, and NDSI remote sensing characteristic curves generated for each plot, the maximum and minimum normalized vegetation index MMNDVI, the maximum and minimum normalized building index MMNDBI, and the maximum and minimum normalized soil index MMNDSI of the plot in each phase are obtained respectively by using the maximum and minimum normalization; The expression of MMNDVI characteristic index is as follows: In the formula, MMNDVI represents the maximum and minimum normalized vegetation index, NDVI represents the normalized vegetation index at a certain moment; max1 represents the maximum value in the NDVI time series data, and min1 represents the minimum value in the NDVI time series data; through the maximum and minimum normalization processing, the values ​​of all calculated MMNDVI remote sensing characteristic indexes are between 0 and 1; The expression of MMNDBI characteristic index is as follows: In the formula, MMNDBI represents the maximum and minimum normalized building index, NDBI represents the normalized building index at a certain moment; max2 represents the maximum value in the NDBI time series data, and min2 represents the minimum value in the NDBI time series data; through the maximum and minimum normalization processing, the values ​​of all calculated MMNDBI remote sensing characteristic indexes are between 0 and 1; The MMNDSI characteristic index expression is as follows: In the formula, MMNDSI represents the maximum and minimum normalized soil index, NDSI represents the normalized soil index at a certain moment, max3 represents the maximum value in the NDSI time series data, and min3 represents the minimum value in the NDSI time series data. Through the maximum and minimum normalization processing, the values ​​of all calculated MMNDSI remote sensing characteristic indexes are between 0 and 1. According to the determined temporary land vector, the usage cycle of temporary land and the calculated NDVI, NDBI, and NDSI remote sensing characteristic indexes of each patch, the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of the temporary land category were drawn respectively through mean processing and maximum and minimum normalization processing, and the typical time series sample characteristics of temporary land were obtained; 30 sample points of non-temporary land were randomly selected to draw the MMNDVI, MMNDBI, and MMNDSI time series characteristic curves of non-temporary land throughout the year.

6. The temporary land use interpretation and reclamation verification method according to claim 1 is characterized in that: In step S2, judging whether the land parcel is temporary land based on the DTW distance criterion includes: Calculate the three remote sensing characteristic indexes including NDVI, NDBI and NDSI of the plot to be judged in each phase, and generate the time series characteristic curve of the three remote sensing characteristic indexes including MMNDVI, MMNDBI and MMNDSI of the land patch to be judged; By comparing the DTW distance similarity between the time series characteristic curves of the three remote sensing characteristic indexes, including MMNDVI, MMNDBI and MMNDSI, of the plot to be judged and the time series characteristic curves of typical temporary land and non-temporary land, it is determined whether the plot to be judged is temporary land.

7. The temporary land use interpretation and reclamation verification method according to claim 6 is characterized in that: After obtaining the typical time series sample characteristics of non-temporary land from the time series characteristic curves of the three remote sensing characteristic indexes, the dynamic time warping (DTW) method is used to calculate the similarity of the time series characteristic curves of the three remote sensing characteristic indexes to determine whether the undetermined land parcel is temporary land. For two time series X n ={c1, c2…c n } and Y m ={q1,q2…q m }, the lengths are |X| and |Y| respectively, calculate the similarity of the two time series, and determine the regular path W=W1,W2…W K ; Among them, Max(|X|,|Y|)≤K≤|X|+|Y|, W K The form of is (i, j), i represents the i coordinate in X, j represents the j coordinate in Y; the normalization path W starts from W1=(1,1) and ends at W1=(|X|,|Y|) to ensure that every coordinate in X and Y appears in W; where c n For the first time series X n The value of the nth phase in m For the second time series Y m The value of the mth phase in K is the Kth point in the path, Max(|X|,|Y|) is the larger value of the two time series lengths, and |X|+|Y| is the sum of the two time series lengths; In the regularized path W, i and j of W(i, j) are monotonically increasing, expressed as: W K =(i,j) W K+1 =(i′,j′) i≤i≤i+1, j≤j≤j+1 Where i and j represent the position index of the Kth point in the first time series and the second time series respectively; W K+1 is the K+1th point in the path, i′ and j′ represent the position index of the K+1th point in the first time series and the second time series respectively, i≤i′≤i+1, j≤j′≤j+1 means that the movement of each step in the path is monotonically increasing; The shortest optimal normalization path is the DTW distance, which is: D(i,j)=Dist(i,j)+min[D(i-1,j),D(i,j-1),D(i-1,j-1)] Where D(i, j) is the accumulated DTW distance between the time series at index i and j, Dist(i, j) is the local Euclidean distance between the i-th point in the first time series and the j-th point in the second time series, min[] represents the smallest one among the three possible paths, that is, the minimum cumulative distance from the top, left or diagonal direction, D(i-1, j) is the accumulated DTW distance between the time series at index i-1 and j, D(i, j-1) is the accumulated DTW distance between the time series at index i and j-1, and D(i-1, j-1) is the accumulated DTW distance between the time series at index i-1 and j-1; The DTW distances between the time series characteristic curves of the three remote sensing characteristic indexes of the pending land parcel and the time series characteristic curves of the three remote sensing characteristic indexes of the temporary land and the non-temporary land are calculated respectively. When the DTW distance with the temporary land is smaller, the pending land parcel is judged as temporary land, otherwise it is judged as non-temporary land.

8. The temporary land use interpretation and reclamation verification method according to claim 1 is characterized in that: In step S3, it is determined whether the temporary land has been reclaimed, including: According to the time series characteristic curves of three remote sensing characteristic indices including MMNDVI, MMNDBI and MMNDSI of temporary land, the slopes of the three remote sensing characteristic indices at the expiration of the temporary land use period and the expiration of the reclamation period are calculated respectively; among which, the expiration of the reclamation period refers to the expiration of one year of the use period; the difference in the slopes of the three remote sensing characteristic indices is combined to determine whether the temporary land with the expiration of the reclamation period has completed land reclamation.

9. The temporary land use interpretation and reclamation verification method according to claim 8, characterized in that: The relevant information for judging whether the temporary land with expiration of reclamation period has completed land reclamation includes: for the patches of temporary land with expiration of reclamation period, the three characteristic curve slopes are calculated by using the three remote sensing characteristic indexes at the expiration time of the time series characteristic curve and the expiration time of reclamation period, respectively. The calculation expression is as follows: In the formula, K MMNDVI , K MMNDBI , K MMNDSI are the curvatures of the time series characteristic curves of the remote sensing characteristic indexes MMNDVI, MMNDBI, and MMNDSI; MMNDVI n ,MMNDBI n ,MMNDSI n are the values ​​of MMNDVI, MMNDBI, and MMNDSI remote sensing characteristic indexes on the time series characteristic curve at the time of reclamation expiration. m ,MMNDBI m ,MMNDSI m are the values ​​of the time of expiration of the use period on the time series characteristic curves of the remote sensing characteristic indexes of MMNDVI, MMNDBI and MMNDSI, respectively; n and m are the time of expiration of the reclamation period and the time of expiration of the use period, respectively; The slope values ​​of the time series characteristic curves of three remote sensing characteristic indices, including MMNDVI, MMNDBI and MMNDSI, are combined to determine whether land reclamation is completed; When K is satisfied MMNDVI >δ1,K MMNDBI <δ2,K MMNDSI <δ3 indicates that the land has been reclaimed and the reclamation effect is good; When K is satisfied MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI ≥δ3 indicates that the land reclamation has not been completed on the plot; When K cannot be satisfied at the same time MMNDVI >δ1,K MMNDBI <δ2,K MMNDSI <δ3 or K MMNDVI ≤δ1,K MMNDBI ≥δ2,K MMNDSI When ≥δ3, manual interpretation is performed; Among them, δ1, δ2, and δ3 are the thresholds of the slopes of the three remote sensing characteristic indexes used to determine whether land reclamation has been completed, which are calculated based on 30 temporary land samples that have completed land reclamation. In practical applications, they are obtained by statistical calculation based on temporary land samples in different regions.

10. A temporary land use interpretation and reclamation verification system, characterized in that: The system is used to control the temporary land use interpretation and reclamation verification method according to any one of claims 1 to 9, and the system comprises: The time series remote sensing index sample feature generation module is used to calculate and generate sample features of typical time series characteristic curves of remote sensing characteristic index of temporary land and non-temporary land by using GEE cloud computing platform and Sentinel-2 satellite remote sensing images; The temporary land use judgment module is used to calculate the time series characteristic curve of the remote sensing characteristic index of the land to be judged, and after linear interpolation or linear compression, calculate the DTW similarity distance between the time series characteristic curve of the remote sensing characteristic index of the land to be judged and the typical time series characteristic curves of temporary land and non-temporary land, and judge whether the land is temporary land based on the DTW distance criterion; The land reclamation judgment module is used to calculate the slope of the time series characteristic curve for temporary land that has expired, based on the remote sensing characteristic index at the time of expiration of the use period and the time of expiration of the reclamation period, to determine whether the temporary land has completed land reclamation.

Images