Clear sandy coastline extraction optimization method based on water body index model

Abstract

The invention provides a clear sandy coastline extraction optimization method based on a water body index model. The method comprises the following steps: acquiring multispectral remote sensing image data at a coast in an area to be subjected to coastline extraction; preprocessing the multispectral remote sensing image data; respectively calculating a normalized difference water body index and an improved normalized difference water body index, and generating corresponding grayscale images; adopting a maximum between-class variance method to determine optimal thresholds of the gray level image generated by the normalized difference water body index and the gray level image generated by the improved normalized difference water body index, generating a binary image of each gray level image based on the optimal thresholds, and obtaining a preliminary water body mask image; introducing a mathematical morphology method to optimize the preliminary water body mask image to obtain an optimized water body mask image; and performing smooth reconstruction on the optimized water body mask image by adopting a sliding window polynomial fitting method to obtain a clear sandy coastline image.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of remote sensing information processing and coastal monitoring, and particularly relates to a clear sandy coastline extraction optimization method based on a water body index model. BACKGROUND

[0002] The position change of the coastline plays an important role in sea level rise, coastline evolution and coastal area planning. As the core boundary of the sea and land, the coastline not only reflects the evolution process of the intertidal zone ecosystem, but also affects the allocation of coastal resources and ecological security. Accurate and real-time sandy coastline information is the basic data for spatial planning of the coastal zone, erosion risk assessment, ecological red line demarcation and disaster prevention and mitigation decision-making. Especially under the background of climate warming, the coastline retreats intensively, and high-precision remote sensing methods are urgently needed to provide support. In the sandy, clear and low turbidity coastal environment, the spectral difference between water and land is small, especially during low tide, wet sand, foam and spray are often mistaken for water, resulting in the landward shift of the coastline. The existing water body index (such as NDWI, MNDWI and AWEI) has a significant decrease in accuracy in such an environment, and needs to be optimized accordingly.

[0003] The existing technology mainly has the following problems:

[0004] Fuzzy boundary between water and sandy: In the clear sandy coast, the spectral reflectance difference between water and dry sand is weak: in the green and near-infrared bands, the reflectance of wet sand and shallow water is close, causing traditional NDWI and other indexes to often misjudge the upper part of the intertidal zone wet sand as water, or misjudge the sand bar exposed by low tide as land, making the water-sand boundary systemically deviate, and the real coastline position is exaggerated or compressed, making it difficult to accurately demarcate the high tide line or the average sea level line.

[0005] Remote sensing image noise interference: Cloud shadow, spray white hat, floating ships, near-shore foam and water grass in tidal ditches and other noises appear as isolated bright spots or dark spots in the water body index gray scale image; these pseudo-features form a large number of voids and fragmented blocks after binarization, causing the coastline to be broken and the burr to be numerous, which requires a large amount of manual interactive editing, seriously restricting the efficiency and stability of automatic extraction.

[0006] Insufficient geometric accuracy of coastline boundary: Pixel-level raster vectorization produces jagged boundaries, and a single pixel error in low-to-medium resolution images (such as 30m) can cause a 10-30m horizontal deviation; higher resolution images can reduce the pixel size, but still cannot depict the micro-scale concave-convex of the sandy beach (beach ridge, gully), resulting in obvious geometric deviation from the measured RTK section or unmanned aerial orthographic image, which is difficult to meet the needs of centimeter-level coastal engineering and change monitoring.

[0007] The shoreline has significant seasonal differences: vegetation growth, storm sand transport, and tidal changes make the shoreline location swing tens of meters in four seasons: in spring and summer, vegetation spectrum covers sand surface, in autumn and winter, storms push sand to the intertidal zone; single-phase images can only capture the instantaneous water edge, and cannot distinguish between short-term tidal fluctuations and long-term erosion / accretion trends, causing confusion in the concept of "shoreline" and making it difficult to establish comparable time series data sets across years and seasons. SUMMARY

[0008] The present application aims to solve the problems in the prior art. The present application provides a clear sandy shoreline extraction optimization method based on a water body index model, which can effectively distinguish sandy boundaries, remove remote sensing image noise interference, improve shoreline boundary geometric accuracy, and generate a standard shoreline sequence across seasons.

[0009] Technical scheme: A clear sandy shoreline extraction optimization method based on a water body index model, comprising the following steps:

[0010] Step 1: Obtain multispectral remote sensing image data of the coast in the area to be extracted;

[0011] Step 2: Preprocess the multispectral remote sensing image data, which includes: band resampling, band fusion, and coordinate system unification;

[0012] Step 3: Calculate the normalized difference water body index and the improved normalized difference water body index respectively, and generate the corresponding gray-scale images;

[0013] Step 4: Determine the best threshold value of the gray-scale images generated by the normalized difference water body index and the improved normalized difference water body index using the maximum inter-class variance method, and generate binary images of each gray-scale image based on the best threshold value to obtain a preliminary water mask image;

[0014] Step 5: Introduce a mathematical morphology method to optimize the preliminary water mask image to obtain an optimized water mask image;

[0015] Step 6: Smoothly reconstruct the optimized water mask image using a sliding window polynomial fitting method to obtain a clear sandy shoreline image.

[0016] Further, in step 2, the preprocessing of the multispectral remote sensing image data includes the following specific operations:

[0017] S301: Select a specific resolution band as a reference and resample other bands using a bilinear interpolation method;

[0018] S302: Select a typical visible light band, RGB band, and fuse the bands with consistent resolution to obtain a composite image;

[0019] S303: Coordinate system of all composite images is unified based on the position of the coastline extraction area.

[0020] Further, the normalized difference water index is calculated according to the following formula:

[0021]

[0022] In the formula, Green represents reflectivity of the green band, and NIR represents reflectivity of the near-infrared band;

[0023] The improved normalized difference water index is calculated according to the following formula:

[0024]

[0025] In the formula, Green represents reflectivity of the green band, and MIR represents reflectivity of the infrared band.

[0026] Further, in step 4, the optimal threshold of the gray-scale image generated by the normalized difference water index and the gray-scale image generated by the improved normalized difference water index is determined by using the maximum inter-class variance method, and specifically includes:

[0027] The optimal threshold t is determined according to the following formula: *

[0028]

[0029] In the formula, t represents the optimal threshold;

[0030] In the formula, P0(t), P1(t) represent the occurrence probabilities of C0 and C1 respectively, and are represented as:

[0031]

[0032] In the formula, i represents a gray level, t represents a segmentation threshold, L represents a total number of gray levels of an image, p i represents a probability of occurrence of the gray level i, and ω0(t) represents a cumulative probability of pixels belonging to the category C0 at the threshold t;

[0033] In the formula, μ0(t), μ1(t) represent the gray level means of C0 and C1 respectively, and are represented as:

[0034]

[0035] In the formula, μ T represents a total gray level mean of the image, and is represented as:

[0036]

[0037] ​Further, in step 5, the mathematical morphology method is introduced to optimize the preliminary water mask image, and an optimized water mask image is obtained, and the specific operation includes:

[0038] S601: defining an original image A and a structural element B as a set in X2, which is a set of symmetry about the origin, and a fixed 3*3 square structural element is adopted;

[0039] S602: performing an expansion operation according to the following formula:

[0040]

[0041] S603: performing an erosion operation according to the following formula:

[0042]

[0043] S604: performing an erosion operation first and then an expansion operation according to the following formula:

[0044]

[0045] S605: performing an expansion operation first and then an erosion operation according to the following formula:

[0046]

[0047] In the formula, A represents the preliminary water mask image, B represents a predefined structural element used for operation, x represents a pixel offset, and z represents a two-dimensional coordinate of a pixel. represents a symmetric point of the two-dimensional coordinate of the pixel.

[0048] Further, in step 6, the sliding window polynomial fitting method is adopted to perform smooth reconstruction on the optimized water mask image, and a clear sandy coastline image is obtained, and the specific operation includes:

[0049] S701: extracting a water edge line from the optimized water mask image, and defining that the water edge line is composed of a series of two-dimensional coordinate points, and represented as A length-k sliding window is applied to the point set, and k points extracted in each sliding window are continuous.

[0050] S702: taking a point as a window value, x as an independent variable, and y as a dependent variable;

[0051] S703: fitting the point set in each window by using an m-order polynomial, and the fitting model is:

[0052]

[0053] The vectorized form is represented as y=Xa+ε.

[0054] where x i represents the base function of polynomial fitting, a i represents the polynomial coefficient corresponding to the term of x i , and epsilon represents the fitting residual.

[0055] where the observation value column vector is y = [y1, y2,..., y k ] T The designed matrix is Each row is The polynomial coefficient to be fitted is a = [a0, a1,..., a m ] T The least square solution is: a = (X T X) -1 X T y.

[0056] S704: The prediction value of the fitted polynomial at the center point of the window is taken as the reconstruction value.

[0057] S705: The window is slid by one bit, and the steps of S702-S704 are repeated, and finally a smooth fitting curve is formed as a clear sandy coastline image.

[0058] Advantages: Compared with the prior art, the present application has the following advantages:

[0059] (1) The method of the present application can realize high-precision, season-consistent and full-automatic extraction and dynamic monitoring of clear sandy coastline, and provide high-precision time-series coastline data for coastal erosion evaluation, ecological red line demarcation and engineering planning.

[0060] (2) The present application can realize high-resolution automatic extraction of clear sandy coastline in complex environment, small geometric error, high-precision dynamic monitoring without manual intervention, and seasonally comparable by constructing NDWI+MNDWI double exponential complementary model, adaptive Otsu threshold segmentation, multi-scale morphological filtering, sliding window polynomial refinement and tide correction seasonal decomposition time series modeling. BRIEF DESCRIPTION OF DRAWINGS

[0061] Figure 1 A work flow diagram of a clear sandy coastline extraction optimization method based on a water body index model provided for an embodiment of the present application;

[0062] Figure 2 A gray scale image generated for the normalized difference water body index;

[0063] Figure 3 An improved gray scale image generated for the normalized difference water body index;

[0064] Figure 4To automatically extract a clear sandy shoreline image. DETAILED DESCRIPTION

[0065] The technical solutions of the embodiments are further described in combination with the drawings and examples.

[0066] The embodiment of the present application provides a clear sandy shoreline extraction optimization method based on a water body index model, as shown in the formula: Figure 1 The embodiment of the present application provides a clear sandy shoreline extraction optimization method based on a water body index model, as shown in the formula:

[0067] Step 1: Selection of the study area, including: selecting a certain place as the center to construct a rectangular range of a certain range as the study area according to the natural geographical features. The natural geographical features can be understood as including the coastal landform type, the tidal type, the tidal level change, the tidal range, the nearshore water quality, so as to have good conditions for carrying out multi-period water edge line extraction and precision analysis; the embodiment of the present application adopts Yantai City and Haiyang City in Shandong Province, and constructs a 5km*5km rectangular range as the study area with Haiyang Bathing Beach as the center. The study area belongs to a typical sandy shoreline, the tidal type is semi-diurnal tide, the tidal level change is regular, the tidal range is moderate, the nearshore water body is clear, and the study area has good conditions for carrying out multi-period water edge line extraction and precision analysis.

[0068] Step 2: Obtain multispectral remote sensing image data at the sandy coast in the study area;

[0069] Step 3: Preprocess the multispectral remote sensing image data, including band resampling, band fusion, and coordinate system unification; specifically including the following steps:

[0070] S301: Select a specific resolution of the green band (B3) with a resolution of 10m as a reference, and resample all bands by using a bilinear interpolation method;

[0071] S302: Select typical visible light bands RGB bands (B4, B3, B2), and complete band fusion processing in the remote sensing image processing and analysis software ENVI, that is, fuse the bands with consistent resolution, so as to obtain a composite image;

[0072] S303: Based on the position of the study area, project all images to a specific geodetic datum WGS1984 coordinate system, and use a specific projection method UTM to realize the unification of the coordinate system.

[0073] Step 4: Construct a water body index model, calculate the normalized difference water body index and the improved normalized difference water body index respectively, and generate corresponding gray scale images, Figure 2 The normalized difference water body index generated gray scale image is shown in the formula, Figure 3 The improved normalized difference water body index generated gray scale image is shown in the formula.

[0074] The normalized difference water index is calculated according to the following formula:

[0075]

[0076] In the formula, Green represents reflectivity of a green wave band, and NIR represents reflectivity of a near-infrared wave band.

[0077] The improved normalized difference water index is calculated according to the following formula:

[0078]

[0079] In the formula, Green represents reflectivity of a green wave band, and MIR represents reflectivity of a middle-infrared wave band.

[0080] Step 5: The optimal threshold value of the gray-scale image generated by the normalized difference water index and the gray-scale image generated by the improved normalized difference water index is automatically determined by using the maximum inter-class variance method, and a binary image is generated to obtain a preliminary water mask image.

[0081] In this step, the calculation formula of the optimal threshold value based on the maximum inter-class variance method is as follows:

[0082]

[0083] In the formula, t * is the optimal threshold value;

[0084] In the formula, P0(t) and P1(t) represent the occurrence probabilities of C0 and C1 respectively, and are expressed as follows:

[0085]

[0086] In the formula, i represents a gray level, and takes a value ranging from 0 to L-1, t represents a segmentation threshold value, L represents the total number of gray levels of an image, and p i represents the probability of occurrence of the gray level i, that is, the ratio of the number of pixels at the gray level to the total number of pixels, and ω0(t) represents the cumulative probability of the pixels belonging to the category C0 at the threshold value t.

[0087] In the formula, μ0(t) and μ1(t) represent the mean values of C0 and C1 respectively, and are expressed as follows:

[0088]

[0089] In the formula, μ T represents the overall mean value of an image, and is expressed as follows:

[0090]

[0091] Step 6: Introduce mathematical morphology method to optimize the preliminary water mask image, including opening operation, closing operation, erosion and dilation operation, to remove noise, fill in the broken part and improve the continuity of water boundary.

[0092] In this step, the mathematical morphology method optimizes the preliminary water mask, including the following steps:

[0093] S601: Define the original image A and the structural element B as a set of X 2 , which is a symmetric set about the origin. A fixed 3x3 square structural element is used to effectively realize boundary smoothing, patch removal and broken repair while maintaining processing efficiency;

[0094] S602: Perform dilation operation to expand the water area, fill in small holes, connect adjacent but discontinuous water pixels, and the calculation formula is:

[0095]

[0096] S603: Perform erosion operation to remove small noise points and non-water pixels, shrink the region boundary and eliminate edge burrs, and the calculation formula is:

[0097]

[0098] S604: First erode and then dilate to effectively eliminate small objects and smooth the boundary while preserving the original region shape, perform opening operation, and the calculation formula is

[0099]

[0100] S605: First dilate and then erode to fill in small holes in the region, connect disconnected regions, smooth the profile, perform closing operation, and the calculation formula is:

[0101]

[0102] In the formula, A represents the preliminary water mask image, i.e. the binary image obtained after calculating the water index (such as NDWI or MNDWI) and threshold segmentation. B represents the predefined structural element used for operation, i.e. a 3x3 square template. x represents an element, i.e. the offset of a pixel, and z represents the two-dimensional coordinates of the pixel, represents the symmetric point of the two-dimensional coordinates of the pixel.

[0103] Step 7: Use sliding window polynomial fitting method to smooth the optimized water boundary and eliminate jagged fluctuations to generate continuous and more geometrically accurate shoreline curves, obtaining Figure 4 the clear sandy shoreline image shown in the figure.

[0104] In this embodiment, the optimized water body boundary is reconstructed by a sliding window polynomial fitting method, including the following steps:

[0105] S701: Extract the water edge line from the optimized water body mask image, and define the preliminary extracted water edge line as a series of two-dimensional coordinate points, represented as Apply a sliding window of length k to the point set, and extract k consecutive points in each sliding window as

[0106] S702: Take a point as the window value, and take x as the independent variable and y as the dependent variable;

[0107] S703: Fit the point set in each window with an m-order polynomial, and the fitting model is:

[0108]

[0109] The vectorized form is represented as y = Xa + ε;

[0110] In the formula, x i represents the base function of polynomial fitting; a i represents the polynomial coefficient corresponding to the x i term, and ε represents the fitting residual.

[0111] Wherein, the observation value column vector is y = [y1, y2,..., y k ] T , and the designed matrix is Each row is the polynomial coefficient to be fitted, and a = [a0, a1,..., a m ] T The least squares solution is: a = (X T X) - 1 X T y;

[0112] S704: Take the predicted value of the fitted polynomial at the center point of the window as the reconstruction value;

[0113] S705: Slide the window by one position, repeat steps S702-S704, and finally form a smooth fitting curve as the clear sandy coastline image.

[0114] Step 8: Spatially compare the clear sandy coastline image with the artificial visually interpreted coastline, calculate the integrity, accuracy and quality three indexes based on the curve matching principle, and complete the precision verification.

[0115] In this step, the precision of the sandy coastline is detected based on the curve matching principle, including the following steps:

[0116] S801: Obtain a clear sandy coastline image based on a manual visual interpretation method;

[0117] S802: Take the real coastline as a reference coastline, and construct a buffer zone with a radius of 0.5-1 pixels for the clear sandy coastline image obtained in step 7; and set a buffer radius of 5-10 meters on the clear sandy coastline image obtained by the manual visual interpretation method;

[0118] S803: The part of the automatic coastline falling into the buffer zone of the reference coastline is recorded as TP1, and the unmatched part is recorded as FP; in the reverse judgment, the part of the reference coastline falling into the buffer zone of the automatic coastline is recorded as TP2, and the unmatched part is recorded as FN;

[0119] S804: Coastline extraction result completeness:

[0120]

[0121] Extraction result accuracy:

[0122]

[0123] Comprehensive evaluation of accuracy and completeness of coastline extraction quality:

[0124]

[0125] Step 9: Repeat steps 1-7 based on four-season remote sensing image data to construct a sandy coastline time series dataset for analyzing coastline dynamic change trends and seasonal response characteristics.

[0126] The embodiment of the present application provides a clear sandy coastline extraction optimization method based on a water body index model, obtains a multispectral remote sensing image sequence covering the sandy coastline, constructs normalized difference water body index and improved normalized difference water body index double index images based on each band, performs through processing on the double index images, obtains coastline extraction points by using maximum inter-class variance method threshold segmentation and mathematical morphology optimization, smoothes and reconstructs the extraction points by using a sliding window polynomial fitting, generates a corrected sandy coastline, and is beneficial to improve the geometric precision of the coastline and the reliability of seasonal monitoring.

Claims

1. A method for extracting and optimizing clear sandy coastlines based on a water index model, characterized in that: Includes the following steps: Step 1: Acquire multispectral remote sensing image data of the coastline within the area to be extracted; Step 2: Preprocess the multispectral remote sensing image data, including: band resampling, band fusion, and coordinate system 1. Step 3: Calculate the normalized difference water index and the improved normalized difference water index respectively, and generate the corresponding grayscale images; Step 4: Using the Otsu's method, determine the optimal threshold for the grayscale image generated by the Normalized Differential Water Index and the grayscale image generated by the Improved Normalized Differential Water Index. Based on the optimal threshold, generate the binarized image of each grayscale image to obtain the preliminary water mask image. Step 5: Introduce mathematical morphology methods to optimize the initial water mask image to obtain the optimized water mask image; Step 6: Use the sliding window polynomial fitting method to smoothly reconstruct the optimized water mask image to obtain a clear sandy coastline image.

2. The method for extracting and optimizing clear sandy coastlines based on a water index model according to claim 1, characterized in that: Step 2, the preprocessing of the multispectral remote sensing image data, specifically includes the following operations: S301: Select a band with a specific resolution as a reference and use bilinear interpolation to resample other bands; S302: Select typical visible light RGB bands, fuse bands with the same resolution to obtain a composite image; S303: Based on the location of the area to be extracted from the coastline, unify the coordinate system of all synthetic images.

3. The method for extracting and optimizing clear sandy coastlines based on a water index model according to claim 1, characterized in that: The normalized differential water index is calculated using the following formula: In the formula, Green represents the reflectivity of the green band, and NIR represents the reflectivity of the near-infrared band; The improved normalized differential water index is calculated using the following formula: In the formula, Green represents the reflectivity of the green band, and MIR represents the reflectivity of the infrared band.

4. The method for extracting and optimizing clear sandy coastlines based on a water index model according to claim 1, characterized in that: Step 4, which involves using the maximum inter-class variance method to determine the optimal threshold for the grayscale image generated by the normalized difference water index and the improved normalized difference water index, specifically includes: The optimal threshold t is determined according to the following formula. * : Where is the optimal threshold; In the formula, P0(t) and P1(t) represent the occurrence probabilities of classes C0 and C1, respectively, and are expressed as: In the formula, i represents the gray level, t represents the segmentation threshold, L represents the total number of gray levels in the image, and p i ω0(t) represents the probability of gray level i appearing, and ω0(t) represents the cumulative probability of a pixel belonging to category C0 at threshold t. In the formula, μ0(t) and μ1(t) represent the mean gray values ​​of classes C0 and C1, respectively, and are expressed as: In the formula, μ T The overall grayscale mean of the image is represented as:

5. The method for extracting and optimizing clear sandy coastlines based on a water index model according to claim 1, characterized in that: In step 5, the mathematical morphology method is introduced to optimize the preliminary water mask image, resulting in an optimized water mask image. Specific operations include: S601: Define the original image A and structuring element B as X 2 The set in It is a symmetric set about the origin, using a fixed 3×3 square structure element; S602: Perform the expansion operation according to the following formula: S603: Perform the etching operation according to the following formula: S604: Corrosion precedes expansion according to the following formula: S605: Expansion precedes corrosion according to the following formula: In the formula, A represents the initial water mask image, B represents the predefined structuring element used for manipulation; x represents the offset of one pixel, and z represents the two-dimensional coordinates of the pixel. Represents the symmetrical point of the pixel's two-dimensional coordinates.

6. The method for extracting and optimizing clear sandy coastlines based on a water index model according to claim 1, characterized in that: In step 6, the sliding window polynomial fitting method is used to smoothly reconstruct the optimized water mask image to obtain a clear sandy coastline image. Specific operations include: S701: Extract the water edge line from the optimized water mask image. The water edge line is defined as a series of two-dimensional coordinate points, represented as... Apply a sliding window of length k to the point set, and extract k consecutive points within each sliding window. S702: Using a certain point as the window value, with x as the independent variable and y as the dependent variable; S703: Within each window, an m-order polynomial is used to fit the point set. The fitting model is as follows: The vectorized form is y = Xa + ε; In the formula, x i The basis functions for polynomial fitting, a i x represents i The polynomial coefficients corresponding to the terms, where ε represents the fitting residuals; The observation column vector is y = [y1, y2, ..., y]. k ] T The designed matrix is Every action The coefficients of the polynomial to be fitted are a = [a0, a1, ..., a2]. m ] T The least squares solution is: a = (X T X) -1 X T y; S704: Use the predicted value of the fitted polynomial at the center point of the window as the reconstructed value; S705: Slide the window one position and repeat steps S702 to S704 to finally form a smooth fitting curve as an image of a clear sandy coastline.

Images