Coastal intertidal zone water line remote sensing extraction technology based on tide level data

Abstract

The invention relates to the technical field of coastline extraction, and particularly discloses a coastal intertidal zone waterline remote sensing extraction technology based on tide level data, which comprises the following steps: S1, acquiring a remote sensing image based on an earthplorer and a copernicius website and preprocessing the remote sensing image; s2, extracting an instantaneous waterline of the remote sensing image by using an NDWI index and a Canny operator to obtain an extraction result; s3, according to an extraction result, discretizing the obtained waterline, and by establishing a function relationship between a discrete point location and a tide level, obtaining an average large tide line so as to determine a preliminary coastline; and S4, correcting the coastline through visual interpretation to obtain a final coastline, and carrying out precision test. According to the coastline extraction method, the NDWI index and the Canny algorithm are used, visual interpretation and an automatic extraction method are combined, the waterline is discretized to determine the position of the average large tide line, and finally the more accurate and efficient coastline extraction method is obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of coastline extraction, and specifically to a remote sensing extraction technology for the water edge line of the coastal intertidal zone based on tide level data. Background Technique

[0002] The change of the coastline is a manifestation of the evolution of the coastal zone resource environment ecosystem under the influence of climate change, sea level change and land-sea interaction, and is a concentrated reflection of human influence on and change of the natural environment as well as adaptation to and response to climate change. Understanding the coastline state through remote sensing methods, monitoring and quantifying the change of the coastline are the basis and key for the environmental protection and planning of the coastal zone. Satellite remote sensing technology is not restricted by conditions such as weather, surface, sea conditions, and geographical environment, and has the characteristics of being fast, accurate, and continuously long-term, and can analyze the dynamic change of the coastline in real time, providing favorable conditions for carrying out the research on the change of the coastline.

[0003] At present, the extraction methods of the coastline are mainly divided into two types: visual interpretation and automatic extraction.

[0004] The visual interpretation method is the most direct method for extracting the coastline. During the process of the big waves and high tides constantly hitting the coast, corresponding traces will be left on the coast, such as the sand ridges formed on the sandy coast, the growth boundaries of some salt-tolerant plants on the biomass coastline, and the seawater trace lines on the artificial coast. These are all important bases for judging the position of the coastline. By field investigation and the different characteristics reflected by different ground objects on the remote sensing image, the corresponding interpretation marks are established, and then the coastline can be interpreted through these marks to obtain the coastline in the strict sense.

[0005] Generally speaking, the coastline obtained by the visual interpretation method has high accuracy and is persuasive. Therefore, some studies on the change of the coastline will adopt the visual interpretation method for extraction. However, the visual interpretation method needs to combine the remote sensing image with field investigation, and at the same time, it needs to manually interpret and interpret the shoreline, with a huge workload.

[0006] The automatic extraction method is the most efficient and practical method in the current process of extracting the coastline, mainly including: the threshold segmentation method and the edge detection method.

[0007] The threshold segmentation method; on the image, the ocean and the land have good spectral differences. Therefore, based on the idea of binary classification, a suitable segmentation threshold is obtained through human-computer interaction or local adaptive strategy to separate the seawater from the land, and the coastline is extracted by positioning the water-land boundary position. Due to the influence of factors such as the spatial heterogeneity of the land, the fluctuation of the seawater, and the tide, the resulting uncertainty problems will greatly affect the availability of the correct segmentation threshold, and the current threshold segmentation methods do not consider the influence of such uncertainty problems.

[0008] Edge detection method; The coastline is located at the boundary between seawater and land. There are significant differences in the characteristics of the features on both sides within the neighborhood, with obvious edge features. Therefore, coastline extraction can be equivalent to an edge extraction problem. As a result, commonly used edge detection algorithms (such as Prewitt, Robert, Sobel, and Canny operators, etc.) have also been widely applied to the research of coastline extraction. The edge detection method completely realizes the target based on the gray-scale change of the image itself. Therefore, it has relatively high requirements for the quality of the image itself. At the same time, when facing large areas, areas with complex coastline types, and the spatial heterogeneity of the land, such methods often fail to achieve satisfactory results and are inefficient. Summary of the Invention

[0009] The purpose of the present invention is to provide a remote sensing extraction technology for the waterline of the coastal intertidal zone based on tide level data to solve the problems mentioned in the above background technology.

[0010] To achieve the above purpose, the present invention provides the following technical solution: A remote sensing extraction technology for the waterline of the coastal intertidal zone based on tide level data, including the following steps:

[0011] S1: Obtain remote sensing images based on the earthexplorer and copernicus websites and perform preprocessing;

[0012] S2: Use the NDWI index and the Canny operator to extract the instantaneous waterline of the remote sensing image to obtain the extraction result;

[0013] S3: According to the extraction result, discretize the obtained waterline, and by establishing the functional relationship between the discrete points and the tide level, obtain the mean high water line of the spring tide to determine the preliminary coastline;

[0014] S4: Through visual interpretation, correct the coastline to obtain the final coastline and conduct accuracy verification.

[0015] Preferably, in the S1, the waterline is extracted using remote sensing image data, and S11: Acquisition of remote sensing images is carried out;

[0016] The principle includes: Obtaining Landsat images and Sentinel-2 images;

[0017] Preprocess the satellite images, including radiometric calibration, atmospheric correction, creating regions of interest, image mosaicking, and image fusion;

[0018] For Sentinel-2 images, the Gram-Schmidt tool is used to select the 8A band with a spatial resolution of 10m from the multi-spectral bands as the panchromatic band, and it is stacked with 6 bands with a resolution of 20m through Layer Stacking to form a multi-band image, obtaining a multi-spectral band image with a resolution of 10m.

[0019] Aiming at the problem of too low resolution of Landsat images, the NNDiffuse image fusion method is adopted to improve the image resolution by fusing the multi-spectral image and the panchromatic image.

[0020] Preferably, in S1, the NDWI of the preprocessed remote sensing image is calculated, that is, S12: preprocess the remote sensing image.

[0021] The principle includes:

[0022] NDWI = (Green - NIR) / (Green + NIR)

[0023] In the formula: Green and NIR are the green light band and the near-infrared band respectively. The calculation of NDWI mainly uses the band operation tool in ENVI software, substitutes the formula "(float(b1) - float(b2)) / (float(b1) + float(b2))", and selects the green light band and the near-infrared band to complete the calculation.

[0024] Using ENVI software, the fused image obtained in step S1 and the NDWI image are stacked through Layer Stacking, and segmentation and extraction are performed in eCognition.

[0025] Set the NDWI classification accuracy weight to 2, perform segmentation at different scales, and obtain the optimal result.

[0026] Preferably, in S2, after completing step S21, that is, using the NDWI index to perform band calculation on the remote sensing image, the mean value and the deviation (D) results can be used to well select the most suitable coastal boundary. The principle includes:

[0027]

[0028] In the formula: Assume the image size is m×n, the image gray scale change range is (0, 255), M(x, y) and F(x, y) represent the low-resolution multi-spectral image and the fused image respectively.

[0029] Preferably, after S21 is completed, step S22 is carried out, that is, the Canny operator is used to process the image after NDWI calculation. After the NDWI-processed remote sensing image is subjected to Canny calculation, the continuity of the image can be effectively improved, and there will be no undetected edges or pseudo-edge phenomena. The Canny operator is written in IDL language to form an extended tool for ENVI:

[0030] First, Gaussian filtering is applied to remove image noise to obtain a noise-removed image; the gradient of the denoised image is calculated; after obtaining the magnitude and direction of the gradient, the pixel points in the image are traversed to remove all non-edge points and suppress non-maximum values; after suppressing non-maximum values, the points of the double-threshold edge are screened to obtain the retained points that form an image. The principle includes:

[0031]

[0032] In the formula: G(x,y) is a two-dimensional Gaussian filtering function, and the parameter σ is the standard deviation of the Gaussian function. When performing image filtering, σ is a parameter that controls its smoothing degree. When σ is selected as a smaller value, the filter has better positioning accuracy, but the denoising effect on the image is relatively poor. When σ is selected as a larger value, the situation of the filter is exactly the opposite. When using Gaussian filtering to denoise the image, it is necessary to select appropriate Gaussian filtering parameters according to specific application requirements. The formula is:

[0033] I(x,y) = f(x,y) * G(x, y)

[0034] In the formula: f(x,y) and I(x,y) are the original grayscale image and the filtered image respectively. For the image smoothed by Gaussian filtering, the gradient magnitude and direction of the filtered image need to be calculated using the finite difference of the first derivative within a 2×2 neighborhood. The formula is:

[0035]

[0036] In the formula: Px(i,j) and Py(i,j) are the partial derivatives of any pixel point (i,j) in the image in its x direction and y direction. Thus, the gradient magnitude M(i,j) and direction θ(i,j) at the pixel point (i,j) can be calculated. The calculation method is as follows:

[0037]

[0038] Preferably, in S22, when further processing the image after non-maximum suppression, it is necessary to use the method of connecting edges with high and low thresholds;

[0039] The principle includes: setting two thresholds, high (Th) and low (Tl), and dividing edge pixels into three categories based on the high and low thresholds, among which pixels with gradient amplitude greater than Th are marked as strong pixels, pixels with gradient amplitude greater than Th but less than Tl are marked as weak edge points, and pixels with gradient amplitude less than Tl are marked as non-edge points. To determine whether a weak edge point is an edge point, it is necessary to see whether there is an edge point greater than Th in the neighborhood pixels of the pixel point. If there is, it is an edge point, otherwise it is a non-edge point.

[0040] Preferably, in said S3, step S31 needs to be performed, that is, the waterside elevation is calculated by the tidal linear interpolation model, based on the two-station interpolation water level correction method, the tidal process curves are respectively drawn using the actual observation results or forecast values ​​of two adjacent tide stations in the study area, and then equidistant interpolation is performed between the two water level process curves, so as to obtain the water level process curve at any position between the two tide stations, and the tidal level at any point between the two tide stations is obtained, and the corresponding waterside elevation is obtained on the basis of the known length of the waterside and the distance along the coastline between the two tide stations;

[0041] In step S32, the tide level at the time of satellite imaging is obtained through the principle of tidal harmonic analysis. The tide at any point is decomposed into the superposition of many partial tides using the principle of tidal harmonic analysis. The amplitude and phase angle of each partial tide are then calculated. After correction by astronomical factors, the harmonic constant of the partial tide is obtained, thereby more accurately predicting the tide level in a certain sea area.

[0042] Preferably, after S32, step S33 is performed, that is, discretizing the waterside line, calculating the average spring tide high level, and discretizing the waterside line using the DSAS plug-in in ArcGIS software; establishing a functional relationship between the discrete points and the tide level to obtain the average spring tide high level;

[0043] When the slope changes are small or relatively consistent on the tidal flat cross section, the position of the tidal characteristic point on the plane can be calculated by using spatial data analysis methods;

[0044] Linearly interpolate the tide level value of the control station at the time of remote sensing image imaging to each discrete point of the instantaneous water edge line to obtain the tide level of the discrete point pair located on the instantaneous water edge line segmentation line. The point position corresponding to the high tide level can be deduced by the similarity principle of trigonometric functions. Connect the average high tide points of all sections into a line to form the water edge line at the time of the remote sensing deduced average high tide level.

[0045] The principles include:

[0046]

[0047] Among them, h 0 is the average high tide level; h 2is the tidal level height at the waterline segmentation point; ((X 0 , Y O )) is the longitude and latitude of the mean spring high tide point on the same segmentation line; ((X 2 , Y 2 )) is the longitude and latitude of the intersection point of the waterline and the segmentation line; (a1, b1), (a2, b2) are the two endpoints of the segmentation line; a is the beach slope, α 1 , α 2 are the projections of the beach slope angle in the x and y directions.

[0048] Preferably, in the step S33, when calculating the mean spring high tide line based on multi-temporal waterlines, according to the known elevations and longitudes of waterline points, combined with the tidal level characteristic value information on the segmentation line, the position of the tidal level characteristic points on the plane can be calculated by using spatial data analysis methods. Connecting all the tidal level characteristic points in sequence can calculate the mean spring high tide line. The formula is as follows:

[0049]

[0050]

[0051] where h 0 is the mean spring high tide level; h 2 , h 3 are the tidal level heights at the corresponding segmentation points of two waterlines; (X 2 , Y 2 ), (X 3 , Y 3 ) are the longitudes and latitudes of the intersection points of two waterlines and the segmentation line respectively; (X 0 , Y O ) is the longitude and latitude of the mean spring high tide point on the same segmentation line.

[0052] Preferably, in the step S4, it is necessary to perform step S41, relying on the national standard hydrographic survey method, that is, the visual interpretation marks of the coastline specified in GB / T12327-1998;

[0053] After the step S41 is completed, step S42 is performed to verify the accuracy by combining the measured points and the simulated points;

[0054] Perform result verification between the simulated points of the calculated mean spring high tide line and the field measured points. The approximate relationship formula for result verification is as follows:

[0055]

[0056] In the formula, η is the sample mean error; (X i , Y i ) is the coordinate of the simulated point; (x i .yi ) are the coordinates of the measured points; α is the standard deviation; d i is the actual error value; n is the number of sample points.

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

[0058] The present invention relates to a remote sensing extraction technology for the waterline of the coastal intertidal zone based on tide level data. By combining the NDWI index and the Canny operator to extract the waterline, calculating the waterline elevation through a tide level linear interpolation model, obtaining the tide level at the satellite imaging moment using the principle of tidal harmonic analysis, on this basis, by discretizing the waterline, determining the position of the mean high water level of spring tides, taking it as the preliminary coastline, and further correcting the coastline by means of visual interpretation to obtain the final coastline. This method combines the NDWI and Canny methods, improves the extraction accuracy of the waterline, combines the automatic extraction technology and manual visual interpretation, improves the extraction efficiency and accuracy of the coastline, and thus proposes a method for determining the coastline using the position of the mean high water level of spring tides. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 is a flowchart of the remote sensing extraction technology for the waterline of the coastal intertidal zone based on tide level data provided by an embodiment of the present invention;

[0060] Figure 2 is a flowchart of another remote sensing extraction technology for the waterline of the coastal intertidal zone based on tide level data provided by an embodiment of the present invention;

[0061] Figure 3 is a research location map provided by an embodiment of the present invention;

[0062] Figure 4 is a flowchart of the waterline extraction based on NDWI and Canny provided by an embodiment of the present invention;

[0063] Figure 5 is a map of the mean high water level of spring tides obtained by discretizing the waterline provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0064] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0065] As Figures 1 to 5 shown, the remote sensing extraction technology for the waterline of the coastal intertidal zone based on tide level data in this embodiment includes the following steps:

[0066] S1: Obtain remote sensing images based on the earthexplorer and copernicus websites and perform preprocessing;

[0067] S2: Use the NDWI index and Canny operator to extract the instantaneous water edge line of the remote sensing image to obtain the extraction result;

[0068] S3: According to the extraction result, discretize the obtained water edge line, establish a functional relationship between the discrete points and the tide level, obtain the mean high water line of the spring tide, and thus determine the preliminary coastline;

[0069] S4: Through visual interpretation, correct the coastline to obtain the final coastline and perform accuracy verification.

[0070] Specifically, in S1, remote sensing image data is used to extract the water edge line, and S11: Acquisition of remote sensing images; the principle includes: obtaining Landsat images and Sentinel-2 images. Landsat images mainly come from the Earth Explorer website (https: / / earthexplorer.usgs.gov / ), and Sentinel-2 images mainly come from the Sentinel website (https: / / browser.dataspace.copernicus.eu / / );

[0071] Perform preprocessing on the satellite images, including radiometric calibration, atmospheric correction, creating regions of interest, image mosaicking, and image fusion;

[0072] For Sentinel-2 images, use the Gram-Schmidt (Gram-Schmidt Spectral Sharpening) tool to use the 8A band with a spatial resolution of 10m in the multi-spectral band as the panchromatic band, and perform Layer Stacking band superposition processing with 6 bands with a resolution of 20m to form a multi-band image, and obtain a multi-spectral band image with a resolution of 10m;

[0073] Regarding the problem of too low resolution of Landsat images, adopt the NNDiffuse image fusion method to improve the image resolution by fusing the multi-spectral image and the panchromatic image.

[0074] Furthermore, perform NDWI calculation on the preprocessed remote sensing images in S1, that is, S12, which can more clearly distinguish land and water bodies, reduce the interference of clouds and fog on the results, enhance the edge information of the land-water boundary, and improve the accuracy of water edge line extraction: perform preprocessing of remote sensing images; the principle includes:

[0075] NDWI = (Green - NIR) / (Green + NIR)

[0076] In the formula: Green and NIR represent the green light band and the near-infrared band respectively. The calculation of NDWI mainly uses the band operation tool in ENVI software, substitutes the formula "(float(b1)-float(b2)) / (float(b1)+float(b2))", and selects the green light band and the near-infrared band to complete the calculation;

[0077] Using ENVI software, the fused image obtained in step S1 and the NDWI image are superimposed through Layer Stacking, and segmentation and extraction are performed in eCognition;

[0078] Set the NDWI classification accuracy weight to 2, perform segmentation at different scales, and obtain the optimal result.

[0079] Furthermore, in S2, step S21 is completed, that is, the NDWI index is used to perform band calculation on the remote sensing image, and the mean value and the deviation (D) results can be used to well select the most suitable coastal boundary. The principles include:

[0080]

[0081] In the formula: Assume that the image size is m×n, the image gray level change range is (0, 255), M(x, y) and F(x, y) represent the low-resolution multispectral image and the fused image respectively.

[0082] Furthermore, after S21 is completed, step S22 is carried out, that is, the Canny operator is used to process the image after NDWI calculation. The remote sensing image after NDWI processing is then subjected to Canny calculation, which can effectively improve the continuity of the image and will not appear the phenomenon of undetected edges or pseudo-edges. The Canny operator is written in IDL language to form an extended tool for ENVI:

[0083] First, apply Gaussian filtering to remove image noise to obtain the image with noise removed; calculate the gradient of the denoised image; after obtaining the magnitude and direction of the gradient, traverse the pixel points in the image, remove all non-edge points, and suppress non-maximum values; after suppressing non-maximum value points, screen the points of the double-threshold edge to obtain the retained points to form an image. The principles include:

[0084]

[0085] Where: G(x, y) is a two-dimensional Gaussian filtering function, and the parameter σ is the standard deviation of the Gaussian function. When performing image filtering, σ is a parameter that controls its smoothness. When σ takes a smaller value, the filter has better positioning accuracy, but the denoising effect on the image is relatively poor. On the contrary, when σ takes a larger value, the situation of the filter is exactly the opposite. When using Gaussian filtering to denoise an image, it is necessary to select appropriate Gaussian filtering parameters according to specific application requirements. The formula is:

[0086] I(x, y) = f(x, y) * G(x, y)

[0087] Where: f(x, y) and I(x, y) are the original grayscale image and the filtered image respectively. For the image smoothed by Gaussian filtering, it is necessary to use the finite difference of the first derivative in the 2×2 neighborhood to calculate the gradient amplitude and direction of the filtered image. The formula is:

[0088]

[0089] Where: Px(i, j) and Py(i, j) are the partial derivatives of any pixel point (i, j) in the image in its x direction and y direction. Thus, the gradient amplitude M(i, j) and direction θ(i, j) at the pixel point (i, j) can be calculated. The calculation method is as follows:

[0090]

[0091] Furthermore, in S22, when further processing the image after non-maximum suppression, it is necessary to use the method of connecting edges with high and low thresholds;

[0092] The principle includes: setting two thresholds, a high threshold (Th) and a low threshold (Tl). Based on the high and low thresholds, the edge pixels are divided into three categories. Among them, the pixel points with a gradient amplitude greater than Th are marked as strong pixel points, the pixel points with a gradient amplitude greater than Th and less than Tl are marked as weak edge points, and the pixel points with a gradient amplitude less than Tl are marked as non-edge points. To determine whether a weak edge point is an edge point, it depends on whether there is an edge point with a gradient amplitude greater than Th among the neighboring pixel points of this pixel point. If there is, it is an edge point; if not, it is a non-edge point.

[0093] Further, in S3, it is necessary to perform step S31, that is, to calculate the waterline elevation through a tidal level linear interpolation model. Based on the two-station interpolation water level correction method, the actual observation results or predicted values of two adjacent tidal stations in the study area are used to draw the tidal level process curves respectively, and then equidistant interpolation is performed between the two water level process curves to obtain the water level process curve at any position between the two tidal stations, and the tidal level at any point between the two tidal stations is obtained. Based on the known waterline length and the shoreline distance between the two tidal stations, the corresponding waterline elevation is obtained; in step S32, the tidal level at the satellite imaging moment is obtained through the principle of tidal harmonic analysis. By applying the principle of tidal harmonic analysis, the tide at any point is decomposed into the superposition of many constituent tides, and then the amplitude and phase angle of each constituent tide are obtained. After correction by astronomical factors, the harmonic constants of the constituent tide are obtained, so as to more accurately predict the tidal level in a certain sea area.

[0094] Further, after S32, step S33 is performed, that is, the waterline is discretized and the mean spring high tide level is calculated. The waterline is discretized by using the DSAS plug-in (Digital Shoreline Analysis System) in the ArcGIS software; a functional relationship is established between the discrete points and the tidal level to obtain the mean spring high tide level.

[0095] In the case where the slope change is small or relatively consistent in the cross-section of the tidal flat, the position of the tidal level characteristic point on the plane can be calculated by using the spatial data analysis method.

[0096] Interpolate the tidal level value of the control station at the remote sensing image imaging moment linearly to each discrete point of the instantaneous waterline, and the tidal level of the discrete point pair located on the instantaneous waterline segmentation line can be obtained. The position of the point corresponding to the high tide level can be calculated by the principle of trigonometric function similarity. Connect the mean spring high tide points of all cross-sections into a line to form the waterline at the moment of the remotely sensed mean spring high tide level.

[0097] The principle includes:

[0098]

[0099] Among them, h 0 is the mean spring high tide level; h 2 is the tidal level height of the waterline segmentation point; ((X 0 , Y O )) are the longitude and latitude of the mean spring high tide point on the same segmentation line; ((X 2 , Y 2 )) are the longitude and latitude of the intersection point of the waterline and the segmentation line; (a1, b1), (a2, b2) are the two endpoints of the segmentation line; a is the shore beach slope, and α 1 , α 2 are the projections of the shore beach slope angle in the x and y directions.

[0100] Further, in step S33, when calculating the mean spring high tide line based on multi-temporal waterlines, according to the known elevations, longitudes and latitudes of waterline points, combined with the tidal level characteristic value information on the dividing line, the spatial data analysis method can be used to calculate the positions of tidal level characteristic points on the plane. Connecting all tidal level characteristic points in sequence can calculate the mean spring high tide line. The formula is as follows:

[0101]

[0102]

[0103] where h 0 is the mean spring high tide level; h 2 , h 3 are the tidal level heights of the corresponding dividing points of two waterlines; (X 2 , Y 2 ), (X 3 , Y 3 ) are the longitudes and latitudes of the intersections of two waterlines and the dividing line respectively; (X 0 , Y O ) are the longitudes and latitudes of the mean spring high tide level points on the same dividing line.

[0104] Furthermore, in step S4, step S41 needs to be carried out, relying on the national standard hydrographic survey method, that is, the visual interpretation marks of the coastline specified in GB / T12327-1998;

[0105] After step S41 is completed, step S42 is carried out to combine the measured points and the simulated points for accuracy verification;

[0106] The result verification is carried out between the simulated points of the calculated mean spring high tide line and the field measured points. The approximate relationship formula for the result verification is as follows:

[0107]

[0108] In the formula, η is the sample mean error; (X i , Y i ) are the coordinates of the simulated points; (x i . y i ) are the coordinates of the measured points; α is the standard deviation; d i is the actual error value; n is the number of sample points.

[0109] Example 1

[0110] Step S1: Obtain remote sensing images based on the earthexplorer and copernicus websites and perform preprocessing;

[0111] In this step, S11: Acquisition of remote sensing images. According to the requirements, Landsat images and Sentinel-2 images are downloaded through the Earth Explorer website

[0112] (https: / / earthexplorer.usgs.gov / ) and the Sentinel website (https: / / browser.dataspace.copernicus.eu / / ) respectively. Since the downloaded Sentinel-2 image data packets are all single-band, band fusion is required to obtain multi-spectral band images.

[0113] In this step, S12: Preprocessing of remote sensing images. The Radiometric Calibration and FLAASH Atmospheric Correction tools in ENVI software are used to perform radiometric calibration and atmospheric correction on the remote sensing images respectively; to meet the needs of the research, ENVI software is used to crop the study area or mosaic multiple scenes of images; the resolution of the remote sensing images is improved through image fusion.

[0114] In this embodiment, as Figure 3 shown, the area with developed tidal creeks in the southern part of the Tiaozini Wetland in Yancheng is selected as the study area. The study area is located along the coast of Yancheng City, Jiangsu Province, China, and is an important transit station and wintering home range on the East Asian-Australasian bird migration route, with longitude and latitude of 120°53′3″~121°1′20″E, 32°37′17″~32°53′8″N. According to the longitude and latitude of the study area, cropping can be directly performed on the Landsat images, while for the Sentinel-2 images, two scenes of images need to be seamlessly mosaicked first and then cropped. To improve the resolution of the remote sensing images, for the Sentinel-2 images, the Gram-Schmidt (Gram-Schmidt Spectral Sharpening) tool is used to use the 8A band with a spatial resolution of 10m as the panchromatic band from the multi-spectral band and perform Layer Stacking band superposition processing with 6 bands with a resolution of 20m to form a multi-band image and obtain a multi-spectral band image with a resolution of 10m; for the problem of too low resolution of the Landsat images, the NNDiffuse image fusion method is adopted. By fusing the multi-spectral image and the panchromatic image, the resolution of the image is improved.

[0115] Step S2: Use the NDWI index and the Canny operator to extract the instantaneous water edge line of the remote sensing image to obtain the extraction result;

[0116] In this step, S21: Calculate the bands of the remote sensing image using the NDWI index. According to the principle of NDWI, select the green band and the near-infrared band to complete the calculation. Use Layer Stacking to overlay the fused image with the NDWI image, and perform segmentation and extraction in eCognition. Visually select the optimal result. Use quantitative statistical methods to evaluate the visual result. The mean refers to the average value of the pixel grayscale, which reflects the average brightness of the image. If the mean is moderate, the visual effect is good. The deviation (D) refers to the difference between the average grayscale value of the original image and the average grayscale value of the fused image. The deviation reflects the average degree of spectral feature change between the fused image and the original multispectral image. The larger the difference, the greater the spectral distortion.

[0117] In this step, S22: Reprocess the image after NDWI calculation using the Canny operator. Implement the Canny operator using IDL language and insert it into the ENVI software as an extended tool. First, use a Gaussian filter to smooth the image to reduce the noise of the water-land boundary line processed by NDWI. Perform convolution operation on the original grayscale image with the Gaussian function to complete the filtering of the image. Then calculate the gradient magnitude and direction of the image. Next is non-maximum suppression, which retains the maximum value of the gradient intensity at each pixel point and filters out other values. Finally, select the points with double-threshold edges to obtain the retained points that form the image.

[0118] In this embodiment, as Figure 4 , use the band operation tool in the ENVI software for the preprocessed remote sensing image, substitute the formula "(float(b1)float(b2)) / (float(b1)+float(b2))", and perform a grayscale histogram analysis on the obtained NDWI image. The result shows that the grayscale histogram of the NDWI image presents a bimodal distribution characteristic. The part of the small wave peak is the land area, and the part of the large wave peak is the water body. The wave valley point is the water-land separation threshold, and the water-land separation threshold is between [-0.2, 0]. Finally, determine the threshold as -0.062. To improve the classification accuracy of water and non-water, in this paper, the preprocessed image and the NDWI image are loaded into the eCognition8.9 software to participate in multi-scale segmentation. Set the NDWI classification accuracy weight to 2, perform segmentation at different scales, and use quantitative statistical methods to evaluate the visual result. The result shows that the mean is moderate and the deviation (D) is large.

[0119] In this embodiment, as Figure 4, the Canny plugin is used to calculate the grayscale image processed by NDWI. First, the weighted average of the pixels around the pixel point is calculated through the filter to obtain the final filtering result. For the Gaussian filter, the closer the point is to the center, the greater the weight value. Of course, the Gaussian filter (Gaussian kernel) is not fixed, and the size of the filter is also variable. The size of the kernel function of the filter is extremely important for the sensitivity of edge detection. Generally speaking, the larger the kernel function, the less likely the edge detection is affected by edge noise, and the more accurate the edge detection information is. Therefore, an important process of edge detection is to determine the size of the Gaussian filter. However, the size of the kernel function often affects the positioning of edge detection, and the positioning accuracy increases with the increase of the kernel function. In this study, a 5×5 kernel function can meet the needs of edge detection. Secondly, the gradient is calculated. The direction of the gradient is always perpendicular to the edge. Usually, the nearby values are taken as horizontal (left, right), vertical (up, down), diagonal (upper right, lower right, upper left, lower left) and other 8 different directions. Therefore, when calculating the gradient, we will get two values: the magnitude and the angle (representing the direction of the gradient) of the gradient. Then, non-maximum suppression is performed. After obtaining the magnitude and direction of the gradient, traverse the pixel points in the image to remove all non-edge points. In the specific implementation, traverse the pixel points one by one, judge whether the current pixel point is the maximum value with the same gradient direction among the surrounding pixel points, and decide whether to suppress the point according to the judgment result. As can be seen from the above description, this step is the edge refinement process. For each pixel point: if the point is the local maximum value in the positive / negative gradient direction, keep the point; if not, suppress the point (set it to zero). Finally, double-threshold edge screening is performed. After completing Gaussian filtering, gradient calculation, and non-maximum suppression, in addition to strong edges, the obtained edge image may also contain some virtual edges. There are various reasons for the existence of virtual edges. It may be generated after edge detection or caused by noise. The strong edges in the image are already in the currently obtained edge image. However, some virtual edges may also be in the edge image. These virtual edges may be generated by the real image or caused by noise. For the latter, they must be removed. The specific division is as follows:

[0120]

[0121] In this embodiment, as Figure 4 , the image processed by NDWI and Canny is imported into the ArcGIS 10.4 software, and the surface-to-line tool is used to obtain the waterway boundary line, and the redundant parts are deleted to finally obtain the instantaneous water edge line.

[0122] Step S3: According to the extraction result, discretize the obtained water edge line, and establish a functional relationship between the discrete points and the tide level to obtain the mean high water line of the spring tide;

[0123] In this step, S31: Based on the two-station interpolation water level correction method, the tidal level process curves are respectively drawn using the actual observation results or predicted values of two adjacent tide gauge stations in the study area, and then equidistant interpolation is performed between the two water level process curves to obtain the water level process curve at any position between the two tide gauge stations, and the tidal level at any point between the two tide gauge stations is obtained. Based on the known lengths of the water edge line and the distance along the shoreline between the two tide gauge stations, the elevation of the corresponding water edge line is obtained.

[0124] In this step, S32: Using the principle of tidal harmonic analysis, the tide at any point is decomposed into the superposition of many constituent tides, and then the amplitude and phase angle of each constituent tide are obtained. After correction by astronomical factors, the harmonic constants of this constituent tide are obtained, so as to more accurately predict the tidal level in a certain sea area.

[0125] In this step, S33: Apply the DSAS module in the ArcGIS software platform, that is, the Digital Shoreline Analysis System, to evenly segment the water edge line, use the Feature to point tool to discretize the water edge line into multiple corresponding scatter points, extract the longitude and latitude coordinates of each scatter point, then use the tidal level linear interpolation model to calculate the elevation of the water edge line points, and finally obtain the mean high water line of spring tides through a functional relationship, and use the mean high water line of spring tides as the preliminary shoreline.

[0126] In this embodiment, as Figure 5 , based on the two-station interpolation water level correction method, the Dafeng Port tide gauge station and the Jianggang tide gauge station are selected as the observation points, and the tidal level process curves are respectively drawn in combination with the tidal level data of the two stations released by the National Marine Science Center. Several short lines parallel to the depth datum are drawn in the middle of the high and low tides of the two tidal level process lines. On each short line, the short line is divided into K equal parts along the direction of A and B. At the high and low tides of the tidal level curve, the high and low tides of the two tide gauge stations are connected into a straight line. Then, near the high and low tides where the curvature is relatively large, the equal division connection lines between the two curves should be parallel to the high (or low) tide connection line. The subsequent equal division connection lines gradually change from being parallel to the high and low tide connection lines to being parallel to the depth datum. The same K equal division is performed on the connection lines, and the corresponding points on each equal division line are connected into a smooth curve, which is the tidal level curve of the interpolated water level zone. Substitute the length of the water edge line and the distance along the shoreline between the two tide gauge stations to obtain the elevation of the corresponding water edge line.

[0127] In this embodiment, as Figure 5, in the actual harmonic analysis process, there are many astronomical tides, but most of their amplitudes are very small and contribute little to the tide level prediction. On the contrary, too many tidal components and truncation errors may reduce its accuracy. In actual analysis, usually a limited number of tidal components that contribute more to the prediction results are selected, which will improve the accuracy of tidal prediction. Therefore, 7 (M2, S2, N2, K1, M4, O1, M6) main tidal components are selected, and the least squares analysis method is used to conduct harmonic analysis on the tide level data to obtain the tide level at the satellite imaging time of the tide gauge station.

[0128] In this embodiment, as Figure 5 , use the ArcGIS software to generate a buffer zone with a distance of 1000m, draw a baseline roughly parallel to the extracted instantaneous waterline, and use it to generate perpendicular lines later. Use the DSAS plug-in (Digital Shoreline Analysis System) in the ArcGIS software. Referring to existing research and based on the characteristics of the instantaneous waterline in this study area, generate perpendicular lines with an interval of 500m and a length of 8000m from the baseline to the side of the instantaneous waterline, and use the Feature topoint tool to discretize the instantaneous waterline into multiple corresponding discrete points, calculate the longitude and latitude, and then obtain the mean high water line of spring tides according to similar triangles.

[0129] Step S4: Through visual interpretation, correct the coastline to obtain the final coastline and conduct result verification;

[0130] In this step, S41: According to the national standard of China "Hydrological Measurement Method" (GB / T 12327-1998), the coastline is "measured based on the actual marks formed during the mean high tide or high tide; the coastline mapping can be determined according to the color, humidity, hardness, soil and vegetation of the coast, as well as alluvial sediments such as driftwood, aquatic plants, and shells. The tidal range in the estuary area is still drawn based on the mean tide and high water line. The coastline can be divided into artificial shorelines and natural shorelines. Artificial shorelines refer to the boundary between the outer side of the tidal flat and artificial traces, which can be identified by remote sensing. The tide usually only has vertical fluctuations on artificial shorelines, without horizontal advance or retreat, mainly including port wharf lines, dike construction lines, salt embankment lines, highway sea dike lines, and river dike embankments; natural shorelines mainly include bedrock shorelines, sandy shorelines, and silt shorelines of silty mud coasts.

[0131] In this step, S42: Use the measured points of the mean high water line of spring tides in the field terrain survey to verify the calculation results of the mean high water line of spring tides in the same year, and use this as the result of the accuracy verification of the tide level characteristic line.

[0132] In this embodiment, according to the standard, the coastline of the study area is divided into artificial shorelines and natural shorelines. The artificial shorelines include salt-cultivation dike shorelines, port wharf lines, riverbank dike shorelines, construction dike shorelines, and road seawall lines; the only type of natural shoreline in the study area is the silty mud coastlines.

[0133] (1) The salt-cultivation dike shoreline is the coastline formed by the dikes of tidal flat aquaculture or salt pans, usually a low-standard earth dike with a tortuous shoreline. The evaporation pond surrounded by the dike has a water color close to that of seawater with low sediment content, and there is generally bare land outside the dike. In the remote sensing image, the layout of the aquaculture area is regular, in a long strip shape. When the cultivation pond is without water, it is similar to bare land, with a high surface reflectivity, showing gray or off-white, and the storage tank has a darker color. The position of the coastline is determined at the outer edge of the aquaculture area.

[0134] (2) The port wharf line is the coastline formed by the permanent coastal structures along the port area. There are residential areas and factories distributed in the port area, and the tone is mostly gray or off-white, which is significantly different from seawater, vegetation, and tidal flats. There are obvious strip-shaped road outlines in this area, and the roads are intricate and easy to interpret. The forms of wharves include coastal wharves, trestle wharves, approach bridge wharves, barge wharves, etc. The coastline of the homeopathic wharf is determined along the front of the wharf, and the coastline of other wharves is determined at the intersection of the terminal root ends.

[0135] (3) The riverbank dike shoreline refers to the coastline of the riverbank, located at the estuary, and is generally parallel to both banks of the river. The river dike coast is similar to the construction of the dike, but its position is at the estuary where the river flows into the sea, and is parallel to both banks of the river. In the image, it has a higher brightness and is white. The position of the riverbank is determined by both sides of the river, and can be determined on the side of the culvert or bridge closest to the sea.

[0136] (4) The construction dike shoreline is an artificial coastline formed by the water-retaining and shore-protecting structures used for the construction of towns, landscapes, or marine leisure and entertainment places. The dike is generally composed of stones, bricks, and concrete, and the image shows high brightness, white, and long strip-shaped. The embankment is usually a dark gray mudflat, showing the characteristics of the filling or construction area. The position of the construction embankment is determined on the sea side of the artificial building.

[0137] (5) The road seawall line is the coastline built along the sea and accessible to vehicles, with the functions of tide protection and wave prevention. It is usually composed of materials such as coastal soil, stones, bricks, and concrete. In the image, it has a high brightness, is white, long and narrow, with an obvious edge line contour, and there is a small amount of bare beach or saline-alkali vegetation outside the embankment. Compared with the salt-containing embankment and dike shorelines, the coastline of the seawall is mostly straight, with a higher construction level and different functions. The position of the road seawall is similar to the dike structure outside the dike.

[0138] (6) The muddy coast is mainly formed by tidal energy, and the muddy coastline is generally located on the muddy coast. On the muddy coast side close to the land, vegetation generally grows vigorously and appears red or dark red in the false color image. On the side close to the sea, there is almost no vegetation.

[0139] In this embodiment, 50 field measurement points were selected within the study area and verified with two types of images respectively. Among them, the distance deviation between the coastline extracted by Sentinel-2 MSI image and the 50 field measurement points is less than 4.71 m. The distance deviation between the coastline extracted by Landsat-8 OLI image and the 50 field measurement points is less than 28.28 m. The average error of the distance between the measured points of the two types of images and the simulated water edge line is less than the maximum allowable error range, the recognition accuracy of the water edge line is high, and the simulation result is accurate.

[0140] Finally, it should be noted that the above are only the preferred embodiments of the present invention and are not used to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, for those skilled in the art, they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A remote sensing extraction technology for coastal intertidal zone water edge based on tidal data, characterized by: The steps include: S1: Obtain remote sensing images based on EarthExplorer and Copernicus websites and perform preprocessing; S2: Use the NDWI index and Canny operator to extract the instantaneous water edge of the remote sensing image and obtain the extraction result; S3: According to the extraction results, the obtained waterside line is discretized, and the average high tide line is obtained by establishing a functional relationship between the discrete points and the tide level, so as to determine the preliminary coastline; S4: Correct the coastline through visual interpretation, obtain the final coastline, and perform accuracy check.

2. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 1 is characterized by: In the S1, the water edge is extracted using the remote sensing image data, and S11: remote sensing image acquisition is performed; The principles include: obtaining Landsat images and Sentinel-2 images; Preprocess satellite images, including radiometric calibration, atmospheric correction, creation of new areas of interest, image mosaicking and image fusion; For Sentinel-2 images, the Gram-Schmidt tool was used to select the 8A band with a spatial resolution of 10m from the multispectral band as the panchromatic band, and layer stacked it with the six bands with a resolution of 20m to form a multi-band image, thus obtaining a multispectral band image with a resolution of 10m. To address the problem of low resolution of Landsat images, the NNDiffuse image fusion method is used to improve the image resolution by fusing multispectral images with panchromatic images.

3. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 2 is characterized by: In the S1, NDWI calculation is performed on the pre-processed remote sensing image, that is, S12: pre-processing the remote sensing image; The principles include: NDWI=Green-NIR / Green+NIR Where: Green and NIR are green light band and near infrared band respectively. The calculation of NDWI mainly uses the band operation tool in ENVI software. Substitute into the formula "(float(b1)-float(b2)) / (float(b1)+float(b2))", select the green light band and near infrared band to complete the calculation; Using ENVI software, the fused image obtained in step S1 was superimposed with the NDWI image through Layer Stacking, and segmented and extracted in eCognition; The NDWI classification accuracy weight was set to 2, and segmentation was performed at different scales to obtain the optimal result.

4. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 3 is characterized by: In step S2, step S21 is completed, that is, the NDWI index is used to calculate the band of the remote sensing image, and the mean value can be used. The sum deviation (D) results are good for selecting the most appropriate coastal boundary. The principles include: In the formula: Assuming that the image size is m×n, the image grayscale range is (0, 255), M(x, y) and F(x, y) represent the low-resolution multispectral image and the fused image respectively.

5. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 4 is characterized in that : After the completion of S21, the step S22 is performed, that is, the image after NDWI calculation is processed by using the Canny operator, and the remote sensing image after NDWI processing is subjected to Canny calculation again, which can effectively improve the continuity of the image, and there will be no undetectable edge or pseudo-edge phenomenon. The Canny operator is written in IDL language to form an extension tool of ENVI: First, apply Gaussian filtering to remove image noise and obtain a denoised image; calculate the gradient of the denoised image; after obtaining the magnitude and direction of the gradient, traverse the pixel points in the image, remove all non-edge points, and suppress non-maximum values; After suppressing the non-maximum points, the points at the double threshold edge are screened to obtain the retained points to form an image. The principles include: Where: G(x, y) is a two-dimensional Gaussian filter function, and the parameter σ is the standard deviation of the Gaussian function. When performing image filtering, σ is a parameter that controls its smoothness. When σ is selected as a small value, the filter has better positioning accuracy, but the denoising effect on the image is relatively poor. When σ is selected as a large value, the situation of the filter is just the opposite. When using Gaussian filtering to denoise the image, it is necessary to select appropriate Gaussian filtering parameters according to specific application requirements. The formula is: I(x,y)=f(xy)*G(x,y) Where: f(x, y) and I(x, y) are the original grayscale image and the filtered image respectively. After the image is smoothed by Gaussian filtering, the finite difference of the first-order derivative in a 2×2 neighborhood is needed to calculate the gradient amplitude and direction of the filtered image. The formula is: Where: Px(i,j) and Py(i,j) are the partial derivatives of any pixel point (i,j) in the image in its x direction and y direction, from which the gradient amplitude M(i,j) and direction θ(i,j) at the pixel point (i,j) can be calculated. The calculation method is as follows:

6. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 5 is characterized by: In S22, when further processing the image after non-maximum suppression, a method of connecting edges with high and low thresholds is required; The principles include: Set two thresholds, high (Th) and low (Tl), and divide edge pixels into three categories based on the high and low thresholds. Among them, pixels with gradient amplitude greater than Th are marked as strong pixels, pixels with gradient amplitude greater than Th but less than Tl are marked as weak edge points, and pixels with gradient amplitude less than Tl are marked as non-edge points. To determine whether a weak edge point is an edge point, we need to see whether there is an edge point with a value greater than Th in the neighborhood pixels of the pixel. If there is, it is an edge point, otherwise it is a non-edge point.

7. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 1 is characterized by: In the above S3, step S31 needs to be performed, that is, the waterside elevation is calculated by the tidal linear interpolation model, based on the two-station interpolation water level correction method, the tidal process curves are drawn respectively by using the actual observation results or forecast values ​​of two adjacent tide stations in the study area, and then equidistant interpolation is performed between the two water level process curves, so as to obtain the water level process curve at any position between the two tide stations, and the tidal level at any point between the two tide stations is obtained, and the corresponding waterside elevation is obtained on the basis of the known length of the waterside and the distance along the coastline between the two tide stations; In step S32, the tide level at the time of satellite imaging is obtained through the principle of tidal harmonic analysis. The tide at any point is decomposed into the superposition of many partial tides using the principle of tidal harmonic analysis. The amplitude and phase angle of each partial tide are then calculated. After correction by astronomical factors, the harmonic constant of the partial tide is obtained, thereby more accurately predicting the tide level in a certain sea area.

8. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 7 is characterized by: After S32, step S33 is performed, that is, the waterside line is discretized, the average spring tide high level is calculated, and the waterside line is discretized using the DSAS plug-in in the ArcGIS software; a functional relationship is established between the discrete points and the tide level to obtain the average spring tide high level; When the slope changes are small or relatively consistent on the tidal flat cross section, the position of the tidal characteristic point on the plane can be calculated by using spatial data analysis methods; Linearly interpolate the tide level value of the control station at the time of remote sensing image imaging to each discrete point of the instantaneous water edge line to obtain the tide level of the discrete point pair located on the instantaneous water edge line segmentation line. The point position corresponding to the high tide level can be deduced by the similarity principle of trigonometric functions. Connect the average high tide points of all sections into a line to form the water edge line at the time of the remote sensing deduced average high tide level. The principles include: Among them, h0 is the average high tide level; h2 is the tidal height of the water edge dividing point; ((X0,Y O )) is the latitude and longitude of the average high tide point on the same dividing line; ((X2, Y2)) is the latitude and longitude of the intersection of the waterside line and the dividing line; (a1, b1) and (a2, b2) are the two endpoints of the dividing line; a is the bank slope, α1 and α2 are the projections of the bank slope angle in the x and y directions.

9. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 7 is characterized by: In the step S33, when the average spring tide high tide line is calculated based on the multi-temporal waterline, the position of the tide feature point on the plane can be calculated using the spatial data analysis method according to the known elevation and longitude and latitude of the waterline point, combined with the tidal feature value information on the dividing line, and all the tidal feature points are connected in sequence to calculate the average spring tide high tide line. The formula is as follows: Among them, h0 is the average high tide level; h2 and h3 are the tidal heights of the two waterside lines corresponding to the dividing point; (X2, Y2) and (X3, Y3) are the longitude and latitude of the intersection of the two waterside lines and the dividing line respectively; (X0, Y O ) is the latitude and longitude of the average high tide point on the same dividing line.

10. The coastal intertidal zone water edge remote sensing extraction technology based on tidal data according to claim 1 is characterized by: In the step S4, the step S41 needs to be performed, relying on the national standard hydrographic method, i.e. GB / T12327-1998, which specifies the visual interpretation signs of the coastline; After the step S41 is completed, the step S42 is performed to combine the measured points with the simulated points to perform accuracy verification; The results of the simulated points of the calculated average spring tide high tide line were verified with the field measured points. The approximate relationship formula of the result verification is as follows: Where η is the sample average error; (X i ,Y i ) is the coordinate of the simulation point; (x i .y i ) is the measured point coordinate; α is the standard deviation; d i is the actual error value; n is the number of sample points.

Images