Large-scale tidal flat range extraction and elevation inversion method and system based on gridding segmentation

Abstract

The invention provides a large-scale tidal flat range extraction and elevation inversion method and system based on gridding segmentation. The method comprises the following steps: S1, screening and preprocessing remote sensing images in a research area to obtain a remote sensing image set; s2, carrying out water body extraction on the screened remote sensing images based on a K-means + + algorithm to obtain a binary image set; s3, summarizing and calculating binarization results to obtain a pixel-by-pixel submerging frequency result, and extracting a tidal flat range; s4, carrying out the inversion of the tidal flat elevation based on the tidal flat elevation inversion model and the grid segmentation to obtain the tidal level information; and S5, spatial filtering post-processing is carried out on the elevation inversion result, and distortion points are removed. Based on the idea of gridding segmentation, the problem that the precision is limited by the assumption that tidal flats with the same tide submerging frequency have the same elevation only in a limited space range during large-scale tidal flat elevation inversion can be solved, and the precision and the speed are both achieved during tidal flat elevation inversion.

Description

Technical Field

[0001] The present invention relates to the technical field of tidal flat research based on remote sensing images, in particular to a large-scale tidal flat range extraction and elevation inversion method and system based on grid segmentation. Background Art

[0002] Tidal flats play an important role in ecological functions, coastal protection, carbon sequestration, natural resources, etc. However, due to their intermediate ecological niche at the junction of land and sea, they are highly dynamic, making them extremely susceptible to various factors, especially human activities and climate change. Against the backdrop of global climate change and rapid socio-economic development, monitoring and studying the changes in tidal flats is a real need and urgency.

[0003] Existing methods for acquiring tidal flat topography are divided into traditional tidal flat topography measurement, lidar measurement, drone photogrammetry, and tidal flat topography inversion methods based on remote sensing. For large-scale tidal flat topography acquisition, satellite remote sensing is currently widely used. Among them, the method of using time-series images to calculate the flooding frequency and then invert the tidal flat topography can capture the dynamic changes of the tidal flat and obtain pixel-by-pixel topographic results, and has been increasingly used. However, this method is generally based on the assumption that "the same flooding frequency corresponds to the same elevation", and existing studies have shown that this assumption is only valid within a limited spatial range. Simply put, when the study area expands, the tides and topographic characteristics within the area will be different, and pixels with the same flooding frequency may not all have the same elevation, which will inevitably reduce the accuracy of the research results.

[0004] Since large-scale tidal flat topography acquisition is currently in high practical demand, it is inevitable to encounter such a contradiction when performing terrain inversion based on inundation frequency: if terrain inversion is performed directly on a larger area, pixels with different elevations but the same inundation frequency will be assigned the same elevation, affecting the accuracy; and if the study area is divided into very detailed small areas, although the terrain inversion accuracy can be improved, it will lead to a complicated implementation process and low scalability of the method.

[0005] Therefore, the present invention has constructed a set of technical processes with simple procedures and high scalability. Based on the idea of ​​grid segmentation, it performs terrain inversion in batches within grids where tides and terrain characteristics do not change much. While retaining the advantage of the ability to depict the details of tidal flat terrain pixel by pixel based on inversion of submergence frequency, it can also solve the above-mentioned contradictions and take into account the speed and accuracy of terrain inversion. It can serve as an effective idea for promoting tidal flat terrain inversion to the national and even global scale. Summary of the Invention

[0006] In view of this, the purpose of the present invention is to provide a large-scale tidal flat range extraction and elevation inversion method and system based on grid segmentation, which can solve the problem that the accuracy of large-scale tidal flat elevation inversion is limited by the "assumption that tidal flats with the same tidal inundation frequency have the same elevation only holds true within a limited spatial range", and achieve both accuracy and speed when performing tidal flat elevation inversion.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: a large-scale tidal flat range extraction and elevation inversion method based on grid segmentation, comprising the following steps:

[0008] Step S1: Acquire Landsat-8 and Sentinel-2 images within a set time and space range, manually screen and remove images with poor image quality or with a large amount of cloud cover in key areas, and perform geometric correction, radiometric correction, atmospheric correction, and cloud removal preprocessing operations to obtain a remote sensing image collection;

[0009] Step S2: extracting water bodies from the filtered remote sensing images based on the K-means++ algorithm to obtain a binary image set;

[0010] Step S3: Summarize and calculate the binarization results to obtain the pixel-by-pixel flooding frequency results and extract the tidal flat range;

[0011] Step S4: Based on the tidal flat elevation inversion model and grid segmentation, the tidal flat elevation is inverted by integrating the inundation frequency results with the tide level information obtained from the TPXO9-atlas tidal model;

[0012] Step S5: Perform spatial filtering on the elevation inversion results to remove distorted points and improve the accuracy of the tidal flat elevation inversion results.

[0013] In a preferred embodiment, step S2 specifically includes the following steps:

[0014] Step S21: Calculate the water index for the remote sensing image. The water index needs to be selected based on the sediment content of the water body in the study area to enhance the water body characteristics and distinguish water bodies from non-water bodies. For turbid water bodies with high sediment content, the modified normalized difference water index (MNDWI) is selected, while for clear water bodies with low sediment content, the automatic water extraction index (AWEI) is selected. The calculation formulas for MNDWI and AWEI are as follows:

[0015]

[0016] Where: ρ Green , ρ SWIR are the spectral band reflectances of the green band and the shortwave infrared band, respectively;

[0017] Automatic Water Extraction Index (AWEI) without shadowsnoshadow :

[0018]

[0019] Automatic Water Extraction Index (AWEI) considering shadows shadow :

[0020]

[0021] Where: ρ Green , ρ Blue , ρ NIR 、 are the reflectances of the spectral bands of green, blue, near infrared, shortwave infrared 1, and shortwave infrared 2, respectively;

[0022] Step S22: clustering the water index calculation results using the K-means++ algorithm. The K-means++ algorithm selects initial cluster centers by moving cluster centers away from each other, and is an improved method of the K-means algorithm.

[0023] Step S23: Select a set number of permanent water body points within the study area, and use the geographic location discrimination method to divide the clustering results into two categories: water bodies and non-water bodies, based on whether there is an intersection with the selected permanent water body points. Set the water body to 1 and the non-water body to 0 to obtain a binary image set.

[0024] In a preferred embodiment: in step S22, the specific steps of the K-means algorithm are as follows:

[0025] Step S221: randomly defining K initial cluster centers;

[0026] Step S222: Calculate the squared Euclidean distance between each data point and the K cluster centers, and assign the data point to the cluster to which the cluster center with the closest Euclidean distance belongs. The calculation formula is as follows;

[0027]

[0028] Among them, d 2 (i,k) represents the squared Euclidean distance between the i-th data point and the k-th cluster center, x ij is the value of the i-th data point in the j-th dimension, s j (k) is the value of the kth cluster center on the jth dimension, and P is the number of dimensions of the variable;

[0029] Step S223: Calculate the centroid of each cluster as the new cluster center of each cluster;

[0030] Step S224: Continuously iterate step 222 and step 223, and stop the iteration when the cluster center no longer changes.

[0031] In a preferred embodiment, the specific steps of the K-means++ clustering algorithm are as follows:

[0032] (1) Randomly select a data point as the first cluster center;

[0033] (2) Calculate the Euclidean distance square D(x) between each remaining sample and the current cluster center i ) 2 ;

[0034] (3) D(x i ) 2 The probability distribution is that the next cluster center is selected according to Formula 5. The farther the point is, the greater the probability of being selected as the new center.

[0035] , n represents the number of data points;

[0036] (4) Repeat steps (2) and (3) until K cluster centers are selected;

[0037] (5) Use these K cluster centers to initialize the cluster centers and then run the K-means algorithm steps.

[0038] In a preferred embodiment, the specific steps of step S3 are as follows:

[0039] Step S31: Comprehensively integrate the binary image set and calculate the flooding frequency of each pixel using the following formula:

[0040]

[0041] Where: P(x,y) represents the flooding frequency of the pixel at a certain location, i' represents the binary image, m is the number of binary images available for each area, B k' (x,y) represents the binarization result of the k'th image;

[0042] Step S32: After calculating the inundation frequency, an empirical threshold of 0.05-0.95 is generally used to extract the tidal flat range. The threshold also needs to be dynamically adjusted according to the actual situation of the study area.

[0043] Step S33: Use vegetation and shoreline masks to remove non-tidal flat factors in the flooding frequency results to obtain the final tidal flat range and the corresponding pixel-by-pixel flooding frequency, and use the "connectedPixelCount" in GEE to calculate the flooding frequency.

[0044] The function removes pixel clusters that are less than or equal to 200 units.

[0045] In a preferred embodiment, the specific steps in step S4 are as follows:

[0046] Step S41: Divide the inundation frequency map into grids with a side length of 10 km, and then determine the tidal point locations for TPXO9-atlas tidal simulation in each grid. The specific method for confirming the simulation point locations is to calculate the centroid of the coordinates of all valued points in each grid and determine whether this centroid is in the valued area. If not, correct it to the coordinates of the valued point closest to the centroid result.

[0047] Step S42: Obtain the relative low tide level H required for inverting the tidal flat elevation l and the relative high tide level H h Specifically, the TPXO9-atlas is used to simulate the hourly tide level within the research time range corresponding to each grid simulation point, and the maximum value is taken as the relative high tide level H h ; Relative low tide level H l The confirmation requires manual review of all the screened remote sensing images, and the transit time of the remote sensing images that are closer to the extracted tidal flat boundary conditions is used as the tide level simulation time of TPXO9-atlas to obtain the relative low tide level H l Since the transit of images is required to obtain relative tide levels, the inversion area needs to be simply divided into tiles according to the smaller tile distribution of Sentinel-2 images. After the inversion is completed in each area, mosaicking is performed to obtain a complete result.

[0048] Step S43: Based on the tidal flat elevation inversion model, each grid is assigned the corresponding relative high tide level H h Relative low tide level H l The tidal flat elevation is calculated, and finally all the small grid calculation results are mosaicked to obtain the elevation inversion result of the entire area.

[0049] In a preferred embodiment, the specific content of the tidal flat elevation inversion model described in step S43 is as follows:

[0050] DEM(x,y)=[(1-P(x,y))×(H h -H l )+H l ];

[0051] Where: DEM (x, y) represents the elevation value of the pixel at position (x, y); P (x, y) represents the frequency of the pixel at position (x, y) being flooded; H h is the relative high tide level; H l Relatively low tide.

[0052] In a preferred embodiment, in step S5, the elevation inversion result is subjected to spatial filtering and post-processing to remove the distorted points. The specific steps are as follows:

[0053] Perform a statistical neighborhood quantile test on the elevation inversion result, calculate the X% quantile of the local N*N window, replace the value of the central pixel with the X% quantile value of the window, and traverse the entire image for filtering.

[0054] The present invention also provides a large-scale tidal flat range extraction and elevation inversion system based on grid segmentation, characterized by running the large-scale tidal flat range extraction and elevation inversion method based on grid segmentation as described in any one of claims 1 to 8.

[0055] Compared with existing technologies, this method offers the following advantages: Based on a grid-based segmentation approach, it addresses the issue of large-scale tidal flat elevation inversion accuracy being limited by the assumption that tidal flats with the same tidal inundation frequency have the same elevation, which only holds true within a limited spatial range. This method achieves both high accuracy and high speed in tidal flat elevation inversion. By leveraging grid-based segmentation, this method improves the accuracy and speed of tidal flat extraction and elevation inversion, providing important data support for coastal research and sustainable development, as well as for decision-making by coastal management agencies.

[0056] The present invention relates to a large-scale tidal flat range extraction and elevation inversion method based on grid segmentation. First, the tidal flat elevation inversion process is effectively solved by the grid segmentation method to solve the problem that the accuracy of large-scale tidal flat elevation inversion is limited by the "assumption that tidal flats with the same tidal inundation frequency have the same elevation only holds true within a limited spatial range". The method combines speed and accuracy. The method is highly operable and easy to implement large-scale tidal flat range extraction and elevation inversion. The method can effectively extract tidal flats in both turbid and clear water areas and has strong scalability. The method is combined with spatial filtering post-processing to eliminate abnormal points, thereby further improving accuracy.

[0057] The present invention clearly aims to solve the problem encountered in large-scale tidal flat elevation inversion, where the accuracy is limited by the assumption that "tidal flats with the same tidal inundation frequency have the same elevation only within a limited spatial range". By using the simple and effective solution of grid segmentation operation, the present invention effectively achieves both high accuracy and speed in large-scale tidal flat elevation inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Attachment Figure 1 It is a schematic flow diagram of the present invention;

[0059] Attachment Figure 2 Schematic diagram of the distribution of images used in the Yangtze River Estuary example and the grid segmentation performed according to the image distribution;

[0060] Attachment Figure 3 This is a detailed schematic diagram of part of the process of the present invention in the Yangtze River Estuary example, taking Chongming Dongtan as an example to illustrate the details;

[0061] Attachment Figure 4 This is the accuracy verification result of the present invention in the Yangtze River Estuary example. DETAILED DESCRIPTION

[0062] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0063] It should be noted that the following detailed descriptions are illustrative and intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present application belongs.

[0064] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application; as used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form, and it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or their combinations.

[0065] Large-scale tidal flat range extraction and elevation inversion method based on grid segmentation, refer to Figure 1-4 , comprising the following steps;

[0066] Step S1: Screening and preprocessing the remote sensing images in the study area to obtain a remote sensing image set;

[0067] Step S2: extracting water bodies from the filtered remote sensing images based on the K-means++ algorithm to obtain a binary image set;

[0068] Step S3: Summarize and calculate the binarization results to obtain the pixel-by-pixel flooding frequency results and extract the tidal flat range;

[0069] Step S4: Based on the tidal flat elevation inversion model and grid segmentation, the tidal flat elevation is inverted by integrating the inundation frequency results with the tide level information obtained from TPXO9-atlas;

[0070] Step S5: Perform spatial filtering on the elevation inversion results to remove distorted points and improve the accuracy of the inversion results.

[0071] In step S2, the remote sensing image is binarized into water and non-water bodies and the submergence frequency result is calculated. This is achieved through water body index calculation, K-means++ clustering, water body extraction, and submergence frequency calculation. Specifically, the following steps are included:

[0072] Step S21: Calculate the water index for the remote sensing image. The water index is selected based on the sediment content of the water in the study area to enhance the water characteristics and distinguish water bodies from non-water bodies. For turbid waters with high sediment content, MNDWI is generally selected, while AWEI is selected for clear waters with low sediment content. The calculation formulas for MNDWI and AWEI are as follows:

[0073] The calculation formulas for MNDWI and AWEI are as follows:

[0074]

[0075] Where: ρ Green , ρ SWIR are the spectral band reflectances of the green band and the shortwave infrared band, respectively.

[0076]

[0077] Where: ρ Green , ρ Blue , ρ NIR , They are the reflectances of the spectral bands of green, blue, near infrared, shortwave infrared 1, and shortwave infrared 2, respectively.

[0078] Step S22: clustering the water index calculation results using the K-means++ algorithm. This algorithm selects initial cluster centers by moving cluster centers away from each other, and is an improved method of the K-means algorithm.

[0079] Step S23: Select a certain number of permanent water body points within the study area. Using the geographic location discrimination method, the clustering results are divided into two categories: water bodies and non-water bodies, based on whether they intersect with the selected permanent water body points. The water body is set to 1 and the non-water body is set to 0, thus obtaining a binary image set.

[0080] In step S22, the specific steps of the K-means algorithm are as follows:

[0081] (1) Randomly define K initial cluster centers;

[0082] (2) Calculate the squared Euclidean distance between each data point and the K cluster centers, and assign the data point to the cluster to which the cluster center with the closest Euclidean distance belongs. The calculation formula is as follows;

[0083]

[0084] Among them, d 2 (i,k) represents the squared Euclidean distance between the i-th data point and the k-th cluster center, x ij is the value of the i-th data point in the j-th dimension, sj (k) is the value of the kth cluster center on the jth dimension, and P is the number of dimensions of the variable. The present invention selects a water body index for clustering, so P is 1.

[0085] (3) Calculate the centroid of each cluster as the new cluster center of each cluster;

[0086] (4) Continuously iterate steps (2) and (3) and stop iterating when the cluster center no longer changes.

[0087] The specific steps of the K-means++ clustering algorithm are as follows:

[0088] (1) Randomly select a data point as the first cluster center;

[0089] (2) Calculate the Euclidean distance square D(x) between each remaining sample and the current cluster center i ) 2 ;

[0090] (3) D(x i ) 2 is the probability distribution. The next cluster center is selected according to Formula 5. The farther the point is, the greater the probability of being selected as the new center.

[0091]

[0092] (4) Repeat steps (2) and (3) until K cluster centers are selected;

[0093] (5) Use these K cluster centers to initialize the cluster centers and then run the K-means algorithm steps.

[0094] In step S3, the method for calculating the pixel-by-pixel flooding frequency result and extracting the tidal flat range is as follows:

[0095] Step S31: Comprehensively integrate the binary image set and calculate the flooding frequency of each pixel using the following formula:

[0096]

[0097] Where: P(x,y) represents the flooding frequency of the pixel at a certain location, i' represents the binary image, m is the number of binary images available for each area, B k' (x,y) represents the binarization result of the k'th image.

[0098] Step S32: After calculating the inundation frequency, an empirical threshold of 0.05-0.95 is generally used to extract the tidal flat range. The threshold also needs to be dynamically adjusted according to the actual situation of the study area.

[0099] Step S33: Use vegetation and shoreline masks to remove non-tidal flat factors in the flooding frequency results, obtain the final tidal flat range and the corresponding pixel-by-pixel flooding frequency, and use the “connectedPixel

[0100] The "Count" function removes pixel clusters that are less than or equal to 200 units.

[0101] In step S4, the grid-based inversion of tidal flat elevation solves the problem of "the assumption that tidal flats with the same tidal inundation frequency have the same elevation only holds true within a limited spatial range" when inverting large-scale tidal flat elevation. This approach achieves both accuracy and speed when inverting tidal flat elevation. The steps for grid-based inversion of tidal flat elevation are as follows:

[0102] Step S41: Divide the inundation frequency map into grids with a side length of 10 km, and then determine the tidal point locations for TPXO9-atlas tidal simulation in each grid. The specific method for confirming the simulation point locations is to calculate the centroid of the coordinates of all valued points in each grid and determine whether this centroid is in the valued area. If not, correct it to the coordinates of the valued point closest to the centroid result.

[0103] Step S42: Obtain the relative low tide level H required for inverting the tidal flat elevation l and the relative high tide level H h Specifically, the TPXO9-atlas is used to simulate the hourly tide level within the research time range corresponding to each grid simulation point, and the maximum value is taken as the relative high tide level H h The confirmation of the relative low tide level Hl requires manual inspection of all the screened remote sensing images, and the transit time of the remote sensing image that is closest to the extracted tidal flat boundary is used as the tide level simulation time of TPXO9-atlas to obtain the relative low tide level Hl. l Since the transit of images is required to obtain relative tide levels, the inversion area needs to be simply divided into tiles according to the tile distribution of the smaller Sentinel-2 images. After the inversion is completed in each area, mosaicking is performed to obtain the complete result.

[0104] Step S43: Based on the tidal flat elevation inversion model, each grid is assigned the corresponding relative high tide level H h Relative low tide level H l The tidal flat elevation is calculated, and finally all the small grid calculation results are mosaicked to obtain the elevation inversion result of the entire area.

[0105] The specific content of the tidal flat elevation inversion model described in step S43 is as follows:

[0106] DEM(x,y)=[(1-P(x,y))×(H h -Hl )+H l ]Formula 7;

[0107] Where: DEM (x, y) represents the elevation value of the pixel at position (x, y); P (x, y) represents the frequency of the pixel at position (x, y) being flooded; H h is the relative high tide level; H l Relatively low tide.

[0108] In step S5, the elevation inversion result is subjected to spatial filtering and post-processing to remove the distorted points. The specific steps are as follows:

[0109] The elevation inversion results are statistically tested for neighborhood quantiles. The 10% quantile of the local 7*7 window is calculated. The value of the central pixel is replaced by the X% quantile value of the window, and the entire image is filtered.

[0110] Specifically, Landsat-8 and Sentinel-2 images available over the four-year period from January 1, 2017, to December 31, 2020, of the Ajiangkou were used.

[0111] like Figure 2 The left figure shows the distribution of images used in this embodiment and the grid segmentation diagram based on the image distribution. As can be seen from the figure, Landsat-8 uses images with two tile positions, and Sentinel-2 uses images with four tile positions. The Yangtze River Estuary was divided into four regions for operation based on the distribution of Sentinel-2 images with smaller image tiles in the Yangtze River Estuary. After screening, the total number of available images in the four regions is 361, with 87, 68, 81, and 125 images in the four regions respectively. Figure 2 As shown in the right figure, grid segmentation operations are performed in four areas respectively, and the tide level simulation points in each grid are calculated.

[0112] like Figure 3 The figure below shows a detailed schematic diagram of part of the process for extracting tidal flats and retrieving elevations in this embodiment. Chongming Dongtan is used as an example to illustrate the details. As can be seen from the figure, the proposed method can effectively separate water bodies, extract tidal flats, and invert elevations. Furthermore, the extracted tidal flats exhibit good detail. Experimental results demonstrate the effectiveness of the proposed method for tidal flat extraction and elevation inversion.

[0113] like Figure 4 As shown in FIG. 1 , the accuracy verification result of the elevation inversion result after spatial filtering and removing the distorted points in this embodiment is obtained. The measured data used for verification comes from the drone measurement results of Chongming Dongtan. According to the accuracy verification results, it can be seen that the method of the present invention has high accuracy.

[0114] The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation of the present invention solves the problem that the accuracy of large-scale tidal flat elevation inversion is limited by the "assumption that tidal flats with the same tidal inundation frequency have the same elevation only holds true within a limited spatial range" through the idea of ​​grid segmentation. It has high accuracy and speed and can be applied to large-scale tidal flat extraction and elevation inversion.

[0115] The above description is only a preferred embodiment of the present invention. All equivalent changes and modifications made according to the scope of the patent application of the present invention should fall within the scope of the present invention.

Claims

1. A large-scale tidal flat range extraction and elevation inversion method based on grid segmentation, characterized by: The following steps are involved: Step S1: Acquire Landsat-8 and Sentinel-2 images within a set time and space range, manually screen and remove images with poor image quality or with a large amount of cloud cover in key areas, and perform geometric correction, radiometric correction, atmospheric correction, and cloud removal preprocessing operations to obtain a remote sensing image collection; Step S2: extracting water bodies from the filtered remote sensing images based on the K-means++ algorithm to obtain a binary image set; Step S3: Summarize and calculate the binarization results to obtain the pixel-by-pixel flooding frequency results and extract the tidal flat range; Step S4: Based on the tidal flat elevation inversion model and grid segmentation, the tidal flat elevation is inverted by integrating the inundation frequency results with the tide level information obtained from the TPXO9-atlas tidal model; Step S5: Perform spatial filtering on the elevation inversion results to remove distorted points and improve the accuracy of the tidal flat elevation inversion results.

2. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 1 is characterized by: Step S2 specifically includes the following steps: Step S21: Calculate the water index for the remote sensing image. The water index needs to be selected based on the sediment content of the water body in the study area to enhance the water body characteristics and distinguish water bodies from non-water bodies. For turbid water bodies with high sediment content, the modified normalized difference water index (MNDWI) is selected, while for clear water bodies with low sediment content, the automatic water extraction index (AWEI) is selected. The calculation formulas for MNDWI and AWEI are as follows: Where: ρ Green , ρ SWIR are the spectral band reflectances of the green band and the shortwave infrared band, respectively; Automatic Water Extraction Index (AWEI) without shadows noshadow : Automatic Water Extraction Index (AWEI) considering shadows shadow : Where: ρ Green , ρ Blue , ρ NIR 、 are the reflectances of the spectral bands of green, blue, near infrared, shortwave infrared 1, and shortwave infrared 2, respectively; Step S22: clustering the water index calculation results using the K-means++ algorithm. The K-means++ algorithm selects initial cluster centers by moving cluster centers away from each other, and is an improved method of the K-means algorithm. Step S23: Select a set number of permanent water body points within the study area, and use the geographic location discrimination method to divide the clustering results into two categories: water bodies and non-water bodies, based on whether there is an intersection with the selected permanent water body points. Set the water body to 1 and the non-water body to 0 to obtain a binary image set.

3. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 2 is characterized by: In step S22, the specific steps of the K-means algorithm are as follows: Step S221: randomly defining K initial cluster centers; Step S222: Calculate the squared Euclidean distance between each data point and the K cluster centers, and assign the data point to the cluster to which the cluster center with the closest Euclidean distance belongs. The calculation formula is as follows; Among them, d 2 (i,k) represents the squared Euclidean distance between the i-th data point and the k-th cluster center, x ij is the value of the i-th data point in the j-th dimension, s j (k) is the value of the kth cluster center on the jth dimension, and P is the number of dimensions of the variable; Step S223: Calculate the centroid of each cluster as the new cluster center of each cluster; Step S224: Continuously iterate step 222 and step 223, and stop the iteration when the cluster center no longer changes.

4. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 3 is characterized by: The specific steps of the K-means++ clustering algorithm are as follows: (1) Randomly select a data point as the first cluster center; (2) Calculate the Euclidean distance square D(x) between each remaining sample and the current cluster center i ) 2 ; (3) D(x i ) 2 is the probability distribution. The next cluster center is selected according to Formula 5. The farther the point is, the greater the probability of being selected as the new center. n represents the number of data points; (4) Repeat steps (2) and (3) until K cluster centers are selected; (5) Use these K cluster centers to initialize the cluster centers and then run the K-means algorithm steps.

5. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 1 is characterized by: The specific steps of step S3 are as follows: Step S31: Comprehensively integrate the binary image set and calculate the flooding frequency of each pixel using the following formula: Where: P(x,y) represents the flooding frequency of the pixel at a certain location, i' represents the binary image, m is the number of binary images available for each area, B k' (x,y) represents the binarization result of the k'th image; Step S32: After calculating the inundation frequency, an empirical threshold of 0.05-0.95 is generally used to extract the tidal flat range. The threshold also needs to be dynamically adjusted according to the actual situation of the study area. Step S33: Use vegetation and shoreline masks to remove non-tidal flat factors in the flooding frequency results to obtain the final tidal flat range and the corresponding pixel-by-pixel flooding frequency, and use the "connectedPixelCount" in GEE to calculate the total flooding frequency. The function removes pixel clusters that are less than or equal to 200 units.

6. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 1 is characterized by: The specific steps in step S4 are as follows: Step S41: Divide the inundation frequency map into grids with a side length of 10 km, and then determine the tidal point locations for TPXO9-atlas tidal simulation in each grid. The specific method for confirming the simulation point locations is to calculate the centroid of the coordinates of all valued points in each grid and determine whether this centroid is in the valued area. If not, correct it to the coordinates of the valued point closest to the centroid result. Step S42: Obtain the relative low tide level H required for inverting the tidal flat elevation l and the relative high tide level H h Specifically, the TPXO9-atlas is used to simulate the hourly tide level within the research time range corresponding to each grid simulation point, and the maximum value is taken as the relative high tide level H h ; Relative low tide level H l The confirmation requires manual review of all the screened remote sensing images, and the transit time of the remote sensing images that are closer to the extracted tidal flat boundary conditions is used as the tide level simulation time of TPXO9-atlas to obtain the relative low tide level H l Since the transit of images is required to obtain relative tide levels, the inversion area needs to be simply divided into tiles according to the smaller tile distribution of Sentinel-2 images. After the inversion is completed in each area, mosaicking is performed to obtain a complete result. Step S43: Based on the tidal flat elevation inversion model, each grid is assigned the corresponding relative high tide level H h Relative low tide level H l The tidal flat elevation is calculated, and finally all the small grid calculation results are mosaicked to obtain the elevation inversion result of the entire area.

7. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 6 is characterized by: The specific content of the tidal flat elevation inversion model described in step S43 is as follows: DEM(x,y)=[(1-P(x,y))×(H h -H l )+H l ]; Where: DEM (x, y) represents the elevation value of the pixel at position (x, y); P (x, y) represents the frequency of the pixel at position (x, y) being flooded; H h is the relative high tide level; H l Relatively low tide.

8. The large-scale tidal flat range extraction and elevation inversion method based on grid segmentation according to claim 1 is characterized by: In step S5, the elevation inversion result is subjected to spatial filtering and post-processing to remove the distorted points. The specific steps are as follows: Perform a statistical neighborhood quantile test on the elevation inversion result, calculate the X% quantile of the local N*N window, replace the value of the central pixel with the X% quantile value of the window, and traverse the entire image for filtering.

9. A large-scale tidal flat range extraction and elevation inversion system based on grid segmentation, characterized by Run the large-scale tidal flat range extraction and elevation inversion method based on grid segmentation as described in any one of claims 1 to 8 above.

Images