Automatic tidal flat extraction method fused with irregular hexagonal grid

Abstract

The invention discloses a tidal flat automatic extraction method fused with an irregular hexagonal grid. The method comprises the following steps: step 1, constructing a time sequence remote sensing data set; step 2, establishing a precision verification sample library; step 3, screening spectral indexes; 4, synthesizing an extreme tide level image; 5, generating an irregular hexagonal grid; step 6, carrying out adaptive threshold segmentation; step 7, secondary elimination of confused ground features; step 8, tidal flat range extraction and precision verification; according to the method, based on a time sequence remote sensing image set, pixels at all positions are screened through a specific spectral index maximum value, a new extreme tide level image is formed through combination, spatial heterogeneity of a complex coastal zone is adapted through an irregular hexagonal grid and a bottom-to-top hierarchical merging strategy, local self-adaptive threshold segmentation is achieved in cooperation with a local Otsu algorithm, and a new extreme tide level image is obtained. The optimal grid scale is automatically determined by means of an elbow method, and the confused ground features are secondarily eliminated through the time sequence standard deviation, so that the problem of miscarriage of the confused ground features is effectively solved, and high-precision full-automatic extraction of the tidal flat is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of environmental remote sensing, in particular to a tidal flat automatic extraction method fusing irregular hexagonal grids. BACKGROUND

[0002] Tidal flat, as a transitional zone between land and sea, plays an important role in ecological protection, disaster prevention and marine economic development. However, due to the influence of natural factors and human activities, the global tidal flat area is continuously shrinking, and it is urgent to support resource sustainable management with accurate and efficient remote sensing extraction technology. The existing tidal flat extraction methods have the following shortcomings: 1. The existing methods rely on instantaneous water level data, and the temporal spectral characteristics are not fully mined, making it difficult to capture the dynamic changes of tidal flat; 2. Although the machine learning model used in the existing method uses temporal features, it has the problems of high classification cost and excessive dependence on training samples; 3. The extraction method based on MSIC-OA algorithm uses global threshold segmentation, which cannot adapt to the spatial heterogeneity of complex coastal zones, resulting in serious misclassification of salt ponds, suspended silt and other confused ground objects, and the extraction accuracy is limited; 4. The applicability of the algorithm used in the existing method to different spatio-temporal resolution images is insufficient, which is difficult to meet the diversified monitoring needs. SUMMARY

[0003] The purpose of the present application is to provide a tidal flat automatic extraction method fusing irregular hexagonal grids to solve the problems raised in the background.

[0004] To achieve the above purpose, the present application provides the following technical scheme: a tidal flat automatic extraction method fusing irregular hexagonal grids, comprising the following steps: step one, constructing a time-series remote sensing dataset; step two, establishing a precision verification sample library; step three, screening spectral indices; step four, synthesizing extreme tide level images; step five, generating irregular hexagonal grids; step six, adaptive threshold segmentation; step seven, secondary exclusion of confused ground objects; step eight, tidal flat range extraction and precision verification.

[0005] In the above step one, the time-series remote sensing images of the target area throughout the year are obtained, the remote sensing images are preprocessed, and a cloud-free time-series image set is constructed.

[0006] In the above step two, the time-series tidal station observation data of the target area throughout the year are obtained, combined with the instantaneous tide height captured by the satellite overflight time image, the local lowest tide image is selected by comparing with the lowest tide height throughout the year; the tidal flat edge points and water area edge points are randomly selected along the tidal flat edge and water area edge respectively, and the completely tidal flat points and confused ground object points are also selected, which together constitute the precision verification sample library.

[0007] In the above step three, four spectral indices NDVI, EVI, NDWI, mNDWI in the cloud-free time-series image set are calculated, the separation effects of each index on the tidal flat and water area, the tidal flat and vegetation, and the building land are analyzed by histogram analysis, the spectral index with the best separation degree of the tidal flat and water area is screened out and recorded as index A, and the spectral index with the best separation degree of the tidal flat and vegetation and building land is screened out and recorded as index B;

[0008] In the above step four, the maximum spectral index synthesis method is used, index A is taken as a quality evaluation band, based on the cloud-free time-series image set, the statistical characteristics of the same position pixels in the time series are combined, the 5% quantile is used to synthesize the minimum water area image, recorded as MSIC min , and the 95% quantile is used to synthesize the maximum water area image, recorded as MSIC max ;

[0009] In the above step five, a bottom-up hierarchical merging strategy is adopted, the minimum grid is divided in a complex area with serious spectral confusion, the adjacent minimum grids are merged in a pure tidal flat area, a series of grid radii with 1000-10000 m and a step length of 200 m is set, a relationship curve of the grid radius and the average standard deviation of index A is drawn, the radius corresponding to the elbow point of the curve is selected as the optimal grid scale, and an irregular hexagonal grid completely covering the target area is generated, each grid being an independent hexagonal classification unit;

[0010] In the above step six, the independent hexagonal classification unit is taken as the minimum classification unit, the Otsu algorithm is applied to the MSIC min image to automatically match the optimal threshold value and retain the permanent water area range; the Otsu algorithm is applied to the MSIC max image to remove the vegetation and construction land and retain the maximum water area range; the spatial difference operation of the maximum water area and the permanent water area is performed to obtain a candidate area image containing the tidal flat and part of the confused ground objects;

[0011] In the above step seven, the time-series standard deviation of index A in the cloud-free time-series image set is calculated, the average of the standard deviations in each grid is calculated based on the irregular hexagonal grid of step five, an adaptive standard deviation threshold is set, the area with a standard deviation fluctuation less than the threshold is regarded as an area without specific periodic change and is determined as a non-tidal flat, and the confused ground objects in the candidate area image are excluded;

[0012] In the above step eight, the candidate area processing results of all irregular hexagonal grids are spliced to form a preliminary extraction result of the tidal flat in the target area; a confusion matrix is constructed by combining the precision verification sample library established in step two, the extraction precision is verified by calculating the precision evaluation index, and the final tidal flat distribution data is output.

[0013] Preferably, in step one, the target area annual time series remote sensing images can be obtained through the GEE cloud platform, and the remote sensing images include Landsat 8 images and Sentinel-2 images, the spatial resolution of the Landsat 8 images is 30m, and the time resolution is 16d, the spatial resolution of the Sentinel-2 images is 10m, and the time resolution is 5d.

[0014] Preferably, in step one, the preprocessing includes radiation calibration, atmospheric correction and cloud removal processing.

[0015] Preferably, in step two, the observation data can be obtained from the China Maritime Service Network through Python crawler.

[0016] Preferably, in step two, the easily confused ground points include artificial wetlands and suspended sediment water areas.

[0017] Preferably, in step three, the histogram analysis result shows that the index A is mNDWI and the index B is EVI.

[0018] Preferably, in step three, the NDVI calculation formula is as follows:

[0019]

[0020] The EVI calculation formula is as follows:

[0021]

[0022] The NDWI calculation formula is as follows:

[0023]

[0024] The mNDWI calculation formula is as follows:

[0025]

[0026] Wherein, is the reflectivity of the near-infrared band; is the reflectivity of the red band; is the reflectivity of the green band; is the reflectivity of the blue band; is the reflectivity of the short-wave infrared band.

[0027] Preferably, in step six, the Otsu algorithm determines the optimal threshold value by calculating the inter-class variance of the image gray histogram, realizes the automatic segmentation of water area and land, tidal flat and vegetation, and the inter-class variance calculation formula is as follows:

[0028]

[0029] Wherein, It represents the sum of squares between groups, describing the variance structure of the dataset; Indicates the number of categories; This indicates the band values ​​selected for classifying different categories; Indicates the category.

[0030] Preferably, in step seven, the adaptive standard deviation threshold is determined by statistical analysis of the temporal standard deviation of index A within each grid, and is used to distinguish the periodic inundation characteristics of tidal flats from the stable spectral characteristics of confusing land features.

[0031] Preferably, in step eight, the accuracy evaluation indicators include overall accuracy, precision, recall, and Kappa coefficient.

[0032] Compared with existing technologies, the beneficial effects of this invention are as follows: Based on a time-series remote sensing image set, this invention filters pixels at each location by the maximum value of a specific spectral index, and combines them to form a new extreme tidal level image. It adapts to the spatial heterogeneity of complex coastal zones with an irregular hexagonal grid and a bottom-up layered merging strategy. It also achieves local adaptive threshold segmentation by combining the local Otsu algorithm, automatically determines the optimal grid scale by using the elbow method, and then eliminates confused features a second time by using the temporal standard deviation. This effectively solves the problem of misclassification of confused features and realizes high-precision fully automatic extraction of tidal flats. Attached Figure Description

[0033] Figure 1 This is a flowchart of the method of the present invention;

[0034] Figure 2 This is a map showing the location of the target area in the experimental example.

[0035] Figure 3 The following are the spectral index histograms of the target region in the experimental example; where (a) is the NDVI index histogram, (b) is the EVI index histogram, (c) is the NDWI index histogram, and (d) is the mNDWI index histogram.

[0036] Figure 4 A schematic diagram of the generation of an irregular hexagonal grid; where (a) is a schematic diagram of the spatial connectivity of a rectangular grid, (b) is a schematic diagram of the spatial connectivity of a hexagonal grid, and (c) is a schematic diagram of the bottom-up layer merging process;

[0037] Figure 5 A schematic diagram illustrating the determination of the optimal grid radius using the elbow method;

[0038] Figure 6The images show a comparison of the extraction results of the MSIC-OA algorithm and the present invention in typical confused areas; where (a) is the minimum water area image synthesized based on Sentinel-2 in the building-tidal flat confused area, (a1) is the extraction result of the MSIC-OA algorithm in the building-tidal flat confused area, and (a2) is the extraction result of the present invention in the building-tidal flat confused area; (b) is the minimum water area image synthesized based on Sentinel-2 in the salt field-suspended sediment area, (b1) is the extraction result of the MSIC-OA algorithm in the salt field-suspended sediment area, and (b2) is the extraction result of the present invention in the salt field-suspended sediment area.

[0039] Figure 7 This is a spatial distribution map of tidal flats in Bohai Bay in 2023 based on Sentinel-2 and Landsat 8 data according to the present invention; wherein, (a) is a spatial distribution map of tidal flats output by the present invention based on Sentinel-2 image, (a1) is an enlarged view of region a1 in Figure (a), and (a2) is an enlarged view of region a2 in Figure (a); (b) is a spatial distribution map of tidal flats output by the present invention based on Landsat 8 image, (b1) is an enlarged view of region b1 in Figure (b), and (b2) is an enlarged view of region b2 in Figure (b);

[0040] Figure 8 This is a flowchart of the steps of the present invention.

[0041] The geographical locations of the above-mentioned remote sensing images are all within China's borders, ensuring the integrity of my country's territory, and are used solely for scientific research purposes. Detailed Implementation

[0042] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0043] Please see the appendix Figure 1 -Appendix Figure 8 The present invention provides an embodiment of an automatic extraction method for tidal flats by fusing irregular hexagonal grids, comprising the following steps: Step 1, constructing a time-series remote sensing dataset; Step 2, establishing an accuracy verification sample library; Step 3, screening spectral indices; Step 4, synthesizing extreme tidal level images; Step 5, generating an irregular hexagonal grid; Step 6, adaptive threshold segmentation; Step 7, secondary exclusion of obfuscated ground features; Step 8, extraction and accuracy verification of tidal flat extent.

[0044] In step one above, year-round time-series remote sensing images of the target area are acquired through the GEE cloud platform. The remote sensing images are then radiometrically calibrated, atmospherically corrected, and declouded to construct a cloud-free time-series image set. The remote sensing images include Landsat 8 images and Sentinel-2 images. The Landsat 8 images have a spatial resolution of 30m and a temporal resolution of 16d, while the Sentinel-2 images have a spatial resolution of 10m and a temporal resolution of 5d.

[0045] In step two above, the annual time-tide level observation data of the tidal station in the target area is obtained from the China Maritime Service Network using Python web crawler. The instantaneous tide height captured by the satellite transit time summary image is combined with the annual lowest tide height to screen the local lowest tide image. Tidal flat edge points and water edge points are randomly selected along the tidal flat edge and water edge respectively. At the same time, complete tidal flat points and easily confused feature points are selected to form an accuracy verification sample library. Among them, easily confused feature points include artificial wetlands and water areas with suspended sediment.

[0046] In step three above, four spectral indices—NDVI, EVI, NDWI, and mNDWI—are calculated in the cloudless time-series image set. Histogram analysis is used to analyze the separation effect of each indice on tidal flats and water areas, tidal flats and vegetation, and built-up areas. The spectral index with the best separation between tidal flats and water areas is selected and denoted as index A; the spectral index with the best separation between tidal flats and vegetation / built-up areas is selected and denoted as index B. Histogram analysis shows that index A is mNDWI and index B is EVI. The NDVI calculation formula is as follows:

[0047]

[0048] The EVI calculation formula is as follows:

[0049]

[0050] The formula for calculating NDWI is as follows:

[0051]

[0052] The formula for calculating mNDWI is as follows:

[0053]

[0054] in, The reflectivity is in the near-infrared band; The reflectivity is for the red band. The reflectivity is for the green band; It is the reflectivity of the blue band; It is the reflectivity of the shortwave infrared band;

[0055] In step four above, the maximum spectral index synthesis method is used, with index A as the quality assessment band. Based on a cloudless time-series image set, the minimum water area image is synthesized by combining the statistical characteristics of pixels at the same location in the time series and applying the 5th quantile, denoted as MSIC. min The maximum water area image was synthesized using the 95th percentile, denoted as MSIC. max ;

[0056] In step five above, a bottom-up hierarchical merging strategy is adopted. In complex areas with severe spectral confusion, the smallest scale grid is divided, and adjacent smallest grids are merged in pure tidal flat areas. A series of grid radii of 1000-10000m with a step size of 200m is set, and the relationship curve between grid radius and the average standard deviation of exponent A is plotted. The radius corresponding to the elbow inflection point of the curve is selected as the optimal grid scale, and an irregular hexagonal grid that completely covers the target area is generated. Each grid is an independent hexagonal classification unit.

[0057] In step six above, the Otsu algorithm is applied to MSIC using independent hexagonal classification units as the smallest classification unit. min In the imagery, the optimal threshold is automatically matched while preserving the permanent water area; the Otsu algorithm is applied to MSIC. max In the imagery, a mask is constructed using index B to remove vegetation and built-up areas, retaining the maximum water area. Candidate areas containing tidal flats and some obfuscated features are obtained through spatial difference calculations between the maximum and permanent water areas. The Otsu algorithm determines the optimal threshold by calculating the maximum inter-class variance of the image's grayscale histogram, achieving automatic segmentation of water and land, and tidal flats and vegetation. The inter-class variance calculation formula is as follows:

[0058]

[0059] in, It represents the sum of squares between groups, describing the variance structure of the dataset; Indicates the number of categories; This indicates the band values ​​selected for classifying different categories; Indicates category;

[0060] In step seven above, the temporal standard deviation of index A in the cloudless temporal image set is calculated. Based on the irregular hexagonal grid in step five, the mean standard deviation within each grid is statistically analyzed, and an adaptive standard deviation threshold is set. Areas with standard deviation fluctuations less than the threshold are considered as areas without specific periodic changes and are determined to be non-tidal flats, thus excluding confusing features in the candidate area images. The adaptive standard deviation threshold is determined through statistical analysis of the temporal standard deviation of index A within each grid and is used to distinguish the periodic inundation characteristics of tidal flats from the stable spectral characteristics of confusing features.

[0061] In step eight above, the candidate region processing results of all irregular hexagonal grids are spliced ​​together to form the preliminary extraction results of the target area tidal flats; a confusion matrix is ​​constructed by combining the accuracy verification sample library established in step two, the overall precision, accuracy, recall and Kappa coefficient are calculated to verify the extraction accuracy, and the final tidal flat distribution data is output.

[0062] Experimental example:

[0063] An automatic extraction method for tidal flats by integrating irregular hexagonal grids includes the following steps: Step 1, constructing a time-series remote sensing dataset; Step 2, establishing an accuracy verification sample library; Step 3, screening spectral indices; Step 4, synthesizing extreme tidal level images; Step 5, generating an irregular hexagonal grid; Step 6, adaptive threshold segmentation; Step 7, secondary exclusion of obfuscated ground features; Step 8, extraction and accuracy verification of tidal flat extent.

[0064] In step one above, the annual time-series Landsat 8 and Sentinel-2 images of the Bohai Bay area (37°30'~39°30'N, 117°30'~121°00'E) in 2023 were acquired through the GEE (Google Earth Engine) cloud platform. After radiometric calibration and atmospheric correction, the S2Cloudless algorithm was used to remove clouds, and Landsat 8 cloud-free image sets and Sentinel-2 cloud-free image sets were constructed respectively.

[0065] In step two above, Python web crawlers were used to obtain the 2023 annual observation data of four tidal stations in the Bohai Bay—Tanggu, Qikou, Huanghua Port, and Caofeidian—from the China Maritime Service Network. Combined with the instantaneous tide height captured by satellite transit time imagery, the local lowest tide images were selected by comparing them with the annual lowest tide height. 300 tidal flat edge points were selected along the edge of the tidal flat, 200 water body edge points were selected along the edge of the water body, 1000 complete tidal flat points and 1000 easily confused feature points (300 salt fields, 300 artificial wetlands, and 400 suspended sediment water bodies) were selected to form an accuracy verification sample library containing 2500 verification samples.

[0066] In step three above, four spectral indices—NDVI (Normalized Difference Vegetation Index), EVI (Enhanced Vegetation Index), NDWI (Normalized Difference Water Index), and mNDWI (Modified Normalized Difference Water Index)—were calculated in the Landsat 8 cloudless image set and the Sentinel-2 cloudless image set. Histogram analysis verified that mNDWI showed the best separation between tidal flats and water areas, while EVI showed the best separation between vegetation and tidal flats. Therefore, mNDWI was determined to be index A, and EVI was determined to be index B.

[0067] In step four above, the `imageCollection.qualityMosaic()` function of the Maximum Spectral Index Composite (MSIC) method is used. Using index A as the quality assessment band, and based on a cloudless time-series image set, the minimum water area image is synthesized by combining the statistical characteristics of pixels at the same location in the time series and applying the 5th percentile. This minimum water area image is denoted as MSIC. min The maximum water area image was synthesized using the 95th percentile, denoted as MSIC. max Capture the extreme submersion and exposure states of tidal flats;

[0068] In step five above, a bottom-up hierarchical merging strategy is adopted to automatically divide the smallest-scale grid in the salt field concentration area in the south of Bohai Bay and the densely built area in the north, and merge the grid in the pure tidal flat area in the middle; a series of grid radii of 1000-10000m with a step size of 200m is set, and the relationship curve between grid radius and the average standard deviation of exponent A is plotted. The radius corresponding to the elbow inflection point of the curve is selected as the optimal grid scale, generating an irregular hexagonal grid that completely covers Bohai Bay, with each grid being an independent hexagonal classification unit;

[0069] In step six above, the Otsu Algorithm (OA) is applied to MSIC, using the independent hexagonal classification unit as the smallest classification unit. min In the imagery, the optimal threshold is automatically matched while preserving the permanent water area; the Otsu algorithm is applied to MSIC. max In the imagery, a mask is constructed using index B to remove vegetation and built-up land, while retaining the maximum water area. By calculating the spatial difference between the maximum water area and the permanent water area, candidate area images containing tidal flats and some obfuscated features are obtained.

[0070] In step seven above, the temporal standard deviation of index A in the cloudless temporal image set is calculated. Based on the irregular hexagonal grid in step five, the mean standard deviation within each grid is statistically analyzed, and an adaptive standard deviation threshold is set. Areas with standard deviation fluctuations less than the threshold are considered as areas without specific periodic changes and are determined to be non-tidal flats, thus excluding confusing features in the candidate area images. The adaptive standard deviation threshold is determined through statistical analysis of the temporal standard deviation of index A within each grid and is used to distinguish the periodic inundation characteristics of tidal flats from the stable spectral characteristics of confusing features.

[0071] In step eight above, the candidate region processing results of all irregular hexagonal grids are stitched together to form the preliminary extraction results of the tidal flats in the target area; a confusion matrix is ​​constructed by combining the accuracy verification sample library established in step two, and the overall accuracy of the Sentinel-2 image extraction results is calculated to be 93.60% and the Kappa coefficient is 0.86, while the overall accuracy of the Landsat 8 image extraction results is 89.60% and the Kappa coefficient is 0.79; finally, the spatial distribution of the tidal flats in Bohai Bay in 2023 is determined.

[0072] This invention combines tidal data, selecting images of the lowest tide solely for verifying the accuracy of tidal flat extraction results. The algorithm applies a time-series remote sensing image set, specifying a particular spectral index (mNDWI) as the "quality assessment band," and combines the maximum values ​​of pixels at the same location across time series to form a new instantaneous image. That is, the generated new instantaneous image is not actually captured by satellite imagery, but rather reconstructed based on each individual pixel. Therefore, in the synthesized MSIC image, each pixel represents the extreme tidal state at its own location.

[0073] Based on the above, the advantages of this invention are as follows: When used, it achieves local adaptive threshold segmentation through an irregular hexagonal grid, effectively separating tidal flats from easily confused features such as salt pans, buildings, and suspended sediment. Compared to the traditional MSIC-OA algorithm, the overall accuracy is improved by 8.4%, and the Kappa coefficient increases from 0.63 to 0.86. The hexagonal grid used better conforms to the natural morphology of the tidal flats. Combined with a bottom-up layered merging strategy, it can divide small-scale grids in areas with severe spectral confusion and merge grids in pure tidal flat areas, avoiding spatial gaps and adapting to the spatial heterogeneity of complex coastal zones. The optimal grid size is automatically determined through a data-driven elbow method, and the entire process of threshold segmentation and confusion feature removal requires no manual intervention, achieving fully automated tidal flat extraction. The technical framework is built on the GEE cloud platform and supports Landsat. 8. Sentinel-2 and other multi-source remote sensing data can efficiently process large-scale, long-term remote sensing data, and are suitable for monitoring tidal flats with images of different spatiotemporal resolutions. They can meet diverse monitoring needs such as coastal zone resource management, ecological protection and disaster prevention, and provide stable and reliable technical support for dynamic monitoring of tidal flats.

[0074] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

Claims

1. A method for automatically extracting tidal flats using an irregular hexagonal grid, comprising the following steps: Step 1: Construct a time-series remote sensing dataset; Step 2: Establish an accuracy verification sample library; Step 3: Screen spectral indices; Step 4: Synthesize extreme tidal level images; Step 5: Generate an irregular hexagonal grid; Step 6: Adaptive threshold segmentation; Step 7: Secondary exclusion of obfuscated features; Step 8: Extraction and accuracy verification of tidal flat extent; Its key feature is: In step one above, the time-series remote sensing images of the target area throughout the year are acquired, the remote sensing images are preprocessed, and a set of cloudless time-series images is constructed. In step two above, the annual time-tide level observation data of the tidal station in the target area is obtained. The instantaneous tide height captured by the satellite transit time summary image is combined with the annual lowest tide height to screen the local lowest tide image. Tidal flat edge points and water edge points are randomly selected along the edge of the tidal flat and the edge of the water area, and complete tidal flat points and easily confused ground features are selected to form an accuracy verification sample library. In step three above, the four spectral indices NDVI, EVI, NDWI, and mNDWI in the cloudless time-series image set are calculated. The separation effect of each index on tidal flats and water areas, tidal flats and vegetation, and built-up areas is analyzed by histogram. The spectral index with the best separation degree between tidal flats and water areas is selected and denoted as index A. The spectral index with the best separation degree between tidal flats and vegetation and built-up areas is selected and denoted as index B. In step four above, the maximum spectral index synthesis method is used, with index A as the quality assessment band. Based on a cloudless time-series image set, the minimum water area image is synthesized by combining the statistical characteristics of pixels at the same location in the time series and applying the 5th quantile, denoted as MSIC. min The maximum water area image was synthesized using the 95th percentile, denoted as MSIC. max ; In step five above, a bottom-up hierarchical merging strategy is adopted. In complex areas with severe spectral confusion, the smallest scale grid is divided, and adjacent smallest grids are merged in pure tidal flat areas. A series of grid radii of 1000-10000m with a step size of 200m is set, and the relationship curve between grid radius and the average standard deviation of exponent A is plotted. The radius corresponding to the elbow inflection point of the curve is selected as the optimal grid scale, and an irregular hexagonal grid that completely covers the target area is generated. Each grid is an independent hexagonal classification unit. In step six above, the Otsu algorithm is applied to MSIC using independent hexagonal classification units as the smallest classification unit. min In the imagery, the optimal threshold is automatically matched and the permanent water area is preserved; the Otsu algorithm is applied to MSIC. max In the imagery, a mask is constructed using index B to remove vegetation and built-up land, while retaining the maximum water area. By calculating the spatial difference between the maximum water area and the permanent water area, candidate area images containing tidal flats and some obfuscated features are obtained. In step seven above, the temporal standard deviation of index A in the cloudless temporal image set is calculated. Based on the irregular hexagonal grid in step five, the mean standard deviation within each grid is statistically analyzed. An adaptive standard deviation threshold is set, and areas with standard deviation fluctuations less than the threshold are regarded as areas without specific periodic changes and are determined to be non-tidal flats, thus excluding confusing ground features in the candidate area images. In step eight above, the candidate region processing results of all irregular hexagonal grids are spliced ​​together to form the preliminary extraction results of the tidal flats in the target area; a confusion matrix is ​​constructed by combining the accuracy verification sample library established in step two, the accuracy evaluation index is calculated to verify the extraction accuracy, and the final tidal flat distribution data is output.

2. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step one, the annual time-series remote sensing images of the target area can be obtained through the GEE cloud platform. The remote sensing images include Landsat 8 images and Sentinel-2 images. The spatial resolution of the Landsat 8 images is 30m and the temporal resolution is 16d, while the spatial resolution of the Sentinel-2 images is 10m and the temporal resolution is 5d.

3. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step one, the preprocessing includes radiometric calibration, atmospheric correction, and cloud removal.

4. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step two, the observation data can be obtained from the China Maritime Service Network using a Python web crawler.

5. The method for automatically extracting tidal flats by incorporating irregular hexagonal grids according to claim 1, characterized in that: In step two, easily confused features include artificial wetlands and water areas with suspended sediment.

6. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step three, the histogram analysis results show that index A is mNDWI and index B is EVI.

7. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step three, the NDVI calculation formula is as follows: ， The EVI calculation formula is as follows: ， The formula for calculating NDWI is as follows: ， The formula for calculating mNDWI is as follows: ， in, The reflectivity is in the near-infrared band; The reflectivity is for the red band. The reflectivity is for the green band; It is the reflectivity of the blue band; It is the reflectivity of the shortwave infrared band.

8. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step six, the Otsu algorithm determines the optimal threshold by calculating the maximum inter-class variance of the image grayscale histogram, thereby achieving automatic segmentation of water areas and land, and tidal flats and vegetation. The formula for calculating the inter-class variance is as follows: ， in, It represents the sum of squares between groups, describing the variance structure of the dataset; Indicates the number of categories; This indicates the band values ​​selected for classifying different categories; Indicates the category.

9. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step seven, the adaptive standard deviation threshold is determined by statistical analysis of the time-series standard deviation of index A within each grid, and is used to distinguish the periodic inundation characteristics of tidal flats from the stable spectral characteristics of confusing land features.

10. The method for automatically extracting tidal flats from an irregular hexagonal grid according to claim 1, characterized in that: In step eight, the accuracy evaluation metrics include overall accuracy, precision, recall, and Kappa coefficient.

Images