Tidal flat terrain inversion method and system

Abstract

The invention discloses a tidal flat terrain inversion method and system, and the method comprises the following steps: 1, obtaining a multi-source remote sensing image and ICESat-2 satellite altimeter data under the condition of less clouds or no clouds, and carrying out the preprocessing of the multi-source remote sensing image and ICESat-2 satellite altimeter data; 2, classifying the preprocessed multi-source remote sensing images by using a support vector machine algorithm, and extracting submerging frequency information of the tidal flat; 3, analyzing photon topographic features of ICESat-2 satellite altimeter data, and accurately extracting prior elevation information of the tidal flat by adopting a multi-stage denoising fusion method; step 4, enabling the tidal flat submerging frequency to be in one-to-one correspondence with the prior elevation information, constructing a tidal flat terrain inversion universal model, and mapping the submerging frequency into the tidal flat terrain; and 5, performing precision evaluation on the inverted tidal flat terrain by using the high-precision airborne LiDAR data. The relation between the tidal flat elevation and the submerging frequency is established, the method can be used for inverting the interannual tidal flat terrain, and high-frequency monitoring of the tidal flat terrain is achieved.

Description

Technical Field

[0001] The present invention relates to a tidal flat terrain remote sensing monitoring method and system, in particular to a tidal flat terrain inversion method and system. Background Art

[0002] As transitional zones where ocean and land frequently interact, tidal flats are crucial to the socioeconomic development of coastal areas, providing valuable habitats, potential land resources, and a buffer zone against marine disasters. However, due to the dual influence of natural factors such as tides, storm surges, and sea level rise, as well as human activities such as artificial reclamation and artificial conservation, tidal flats around the world are experiencing varying degrees of reduction and erosion. Furthermore, the location, shape, and extent of tidal flats often undergo dramatic changes, exhibiting strong spatiotemporal dynamics, resulting in numerous challenges for existing tidal flat topography inversion methods. Therefore, there is an urgent need to develop fast and accurate methods for acquiring and dynamically updating tidal flat topography to achieve accurate inversion and dynamic monitoring of tidal flat topography, providing a scientific basis for coastal protection and sustainable development.

[0003] Currently, remote sensing technology has been widely used for tidal flat topography inversion due to its advantages such as wide coverage, strong periodicity, and high timeliness. Among them, the waterline method is the most widely used, but it suffers from problems such as cumbersome procedures, low automation, and low elevation calibration accuracy. The inundation frequency method is an improvement on the waterline method. It calculates the inundation frequency from time-series remote sensing imagery and then calculates the mean elevation of pixels with the same inundation frequency to invert the tidal flat topography. However, existing topographic tidal flat inversion techniques based on the inundation frequency method typically use a single data source, capturing the full tidal flat inundation cycle with a step size of 2-4 years. This results in averaging the tidal flat topography over that period, making it impossible to achieve high-frequency monitoring of the tidal flat topography and difficult to visualize subtle changes in topography between different years. Furthermore, existing technologies have not yet established a functional relationship between inundation frequency and tidal flat elevation, lacking a unified and universal mathematical model. Summary of the Invention

[0004] Purpose of the invention: In order to overcome the deficiencies in the prior art, the purpose of the present invention is to provide a high-precision tidal flat topography inversion method. Another purpose of the present invention is to provide a tidal flat topography inversion system with accurate inversion and dynamic monitoring.

[0005] Technical solution: The tidal flat topography inversion method described in the present invention comprises the following steps:

[0006] Step 1: Acquire multi-source remote sensing images and ICESat-2 satellite altimeter data under light or no cloud conditions and perform preprocessing.

[0007] Step 2: Use the support vector machine algorithm to classify the pre-processed multi-source remote sensing images and extract the flooding frequency information of the tidal flat;

[0008] Step 3: By analyzing the photon topographic characteristics of ICESat-2 satellite altimeter data, a multi-level denoising and fusion method is used to accurately extract the prior elevation information of the tidal flat;

[0009] Step 4: The tidal flat inundation frequency is mapped to the prior elevation information one by one, and a general tidal flat topography inversion model is constructed to map the inundation frequency to the tidal flat topography.

[0010] Step 5: Use high-precision airborne LiDAR data to evaluate the accuracy of the inverted tidal flat topography.

[0011] Furthermore, in step one, the multi-source remote sensing images include Landsat-8 / 9, Sentinel-2, HJ-2A / B and GF-1 optical images. The preprocessing of the multi-source remote sensing images includes radiometric calibration, atmospheric correction, geometric precision correction, and unification of spatial resolution and coordinate system. The preprocessing of ICESat-2 satellite altimeter data includes confidence screening and spatial range screening.

[0012] Furthermore, in step 2, the tidal flat inundation frequency is calculated according to the following formula:

[0013] F I =N sea_water / (N tidal_flat +N sea_water )

[0014] Among them, F I is the flooding frequency of a pixel on the tidal flat, N tidal_flat is the frequency of a certain pixel being observed as a tidal flat, N sea_water is the frequency of observing seawater at a certain pixel on the tidal flat.

[0015] Furthermore, in step 3, the elevation threshold for removing seawater photons is calculated according to the following formula:

[0016] T=H msl +2·σ msl

[0017] Among them, T is the elevation threshold, H msl is the average sea level height of photons obtained from the elevation histogram statistics, σ msl is the standard deviation of seawater photon elevation.

[0018] Furthermore, in step 3, a density clustering method based on the adaptive reachable distance of the Otsu method is used to perform preliminary denoising on the tidal flat photons. The calculation formulas for the reachable distance and the Otsu method are as follows:

[0019] RD(o,o')=max(core_dist(o'),dist(o,o'))

[0020]

[0021] Where RD(o,o′) is the reachable distance from point o to point o′, core_dist(o′) is the core distance of point o′, and dist(o,o′) is the Euclidean distance between point o and point o′; σ B 2 (t) is the inter-class variance under threshold t, ω1(t) is the pixel ratio of the foreground class corresponding to threshold t, and ω2(t) is the pixel ratio of the background class corresponding to threshold t, μ1(t) is the mean of the foreground, μ2(t) is the mean of the background, μ T is the overall mean of the point cloud.

[0022] Furthermore, in step 3, the calculation formula for the kernel density in the Gaussian low-pass filter is as follows:

[0023]

[0024] in, is the position h i The kernel density estimate at h i is the elevation value of the i-th photon point, h j is the elevation value of the j-th photon point, n is the total number of photons involved in the estimation, and σ is the standard deviation of the kernel function, which controls the smoothness of the kernel density curve.

[0025] Furthermore, in step 4, a grid is established based on the pixel size of the tidal flat inundation frequency, and the prior elevation is averaged using grid filtering to obtain aligned prior elevation and inundation frequency data.

[0026] Furthermore, in step 4, the general model for tidal flat topography inversion is as follows:

[0027] H l =wf+b

[0028]

[0029] Among them, H l is a linear model, H nl is a nonlinear model, f is the flooding frequency, q is the highest order of the nonlinear model, w, b, u i , u0 is the model coefficient.

[0030] Furthermore, in step 5, the calculation formula for accuracy evaluation is:

[0031]

[0032] Among them, RMSE is the root mean square error, MAE is the mean absolute error, Bias is the bias, R 2 is the coefficient of determination, yest Reverse the terrain elevation value for regression modeling, y obs is the elevation value of the airborne LiDAR data, is the average elevation value of the airborne LiDAR data, and n is the number of samples.

[0033] The tidal flat topography inversion system of the present invention comprises:

[0034] The inundation frequency calculation module is used to pre-process multi-source remote sensing images, combine the land-water segmentation algorithm to map the tidal flats, and obtain the inundation frequency data of the tidal flats;

[0035] The tidal flat elevation extraction module is used to preprocess and denoise the lidar satellite altimeter data to obtain the prior elevation information of the tidal flat;

[0036] The tidal flat topography inversion module is used to construct linear and nonlinear tidal flat topography inversion models based on regression modeling, and to achieve accurate inversion of tidal flat topography by combining tidal flat inundation frequency.

[0037] Beneficial effects: Compared with the prior art, the present invention has the following significant features:

[0038] 1. The interannual tidal flat inundation frequency was extracted based on multi-source remote sensing images, improving the temporal resolution of the inundation frequency. Furthermore, linear and nonlinear quantitative mapping relationships between tidal flat elevation and inundation frequency were established by coupling regression modeling technology. This enabled the constructed tidal flat topography to be used to invert the interannual tidal flat topography, thus achieving high-frequency monitoring of the tidal flat topography.

[0039] 2. It can solve the limitations of traditional tidal flat topography monitoring technology in high-frequency monitoring, provide an efficient and universal method for dynamic updating of tidal flat topography, and thus provide technical support for coastal protection and sustainable development, which has great practical application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 It is a workflow diagram of the present invention;

[0041] Figure 2 It is a time distribution map of multi-source remote sensing images used in the present invention;

[0042] Figure 3 is a spatial distribution diagram of submergence frequency of the present invention;

[0043] Figure 4 It is a tidal flat topography inversion model constructed by the present invention;

[0044] Figure 5 is a spatial distribution diagram of the tidal flat terrain inversion results of the present invention, wherein (a) is the terrain inversion result based on the linear model, and (b) is the terrain inversion result based on the nonlinear model;

[0045] Figure 6 : is a scatter plot of the accuracy evaluation of the tidal flat topography inversion results of the present invention, wherein (a) is the accuracy evaluation result of the linear model based on ICESat-2 data, (b) is the accuracy evaluation result of the nonlinear model based on ICESat-2 data, (c) is the accuracy evaluation result of the linear model based on airborne LiDAR data, and (d) is the accuracy evaluation result of the nonlinear model based on airborne LiDAR data;

[0046] Figure 7 It is a system framework diagram of the present invention. DETAILED DESCRIPTION

[0047] Example 1

[0048] This example takes the radial sand ridges in northern Jiangsu as an example to conduct a precise inversion experiment on the tidal flat topography in this area. Figure 1 As shown in FIG, a tidal flat topography inversion method based on the collaboration of multi-source remote sensing and regression modeling includes the following steps:

[0049] Step 1: Obtain multi-source remote sensing images and satellite altimeter data covering the study period and target area. This example collected 83 images from Landsat-8 / 9, Sentinel-2, HJ-2A / B, and GF-1 in 2024, as well as 9 ICESat-2ATL03 orbital data. Figure 2 As shown in the figure, multi-source remote sensing images were subjected to radiometric calibration, atmospheric correction, and geometric precision correction, with a unified spatial resolution of 30 meters and a coordinate system of WGS_1984_UTM_Zone_51N. ICESat-2 data preprocessing involves selecting photon points with a confidence level of at least 2 based on the confidence level information of the ATL03 product. Spatial range screening is also performed based on the vector boundary of the northern Jiangsu tidal flats to obtain valid ICESat-2 observation data.

[0050] Step 2: Remote sensing calculation of tidal flat inundation frequency. Based on the multi-source remote sensing image of the Subei radial sand ridge group obtained in step 1, the support vector machine algorithm in ArcGIS Pro was used to perform land and water segmentation to obtain a binary image. Then, Arcpy was used to calculate the tidal flat inundation frequency. Figure 3 As shown in Figure 3, the inundation frequency map obtained based on multi-source remote sensing images can characterize the inundation state of the tidal flat and also reflect the local details of the tidal flat such as the tidal gully.

[0051] The frequency of tidal flat inundation is calculated using the following formula:

[0052] F I =N sea_water / (N tidal_flat +N sea_water )

[0053] Among them, F I is the flooding frequency of a pixel on the tidal flat, N tidal_flat is the frequency of a certain pixel being observed as a tidal flat, N sea_water is the frequency of observing seawater at a certain pixel on the tidal flat.

[0054] Step 3: Extracting the prior elevation of the tidal flat. The specific method is to first perform a preliminary screening based on the ICESat-2 ATL03 photon elevation data obtained in step 1 using the photon confidence attribute, and then remove invalid photons based on the spatial range of the tidal flat. The elevation threshold for removing seawater photons is calculated using the following formula:

[0055] T=H msl +2·σ msl

[0056] Among them, T is the elevation threshold, H msl is the average sea level height of photons obtained from the elevation histogram statistics, σ msl is the standard deviation of seawater photon elevation.

[0057] The elevation statistical histogram method is then used to remove sea surface photons. The OPTICS (Ordering Points To Identify the Clustering Structure) density clustering algorithm is then used to construct a reachable distance map. Photon outliers are adaptively identified based on core distance and neighborhood density, effectively removing low-density noise points.

[0058] The density clustering method based on the adaptive reachable distance of the Otsu method is used to perform preliminary denoising on the tidal flat photons. The calculation formulas for the reachable distance and the Otsu method are as follows:

[0059] RD(o,o')=max(core_dist(o'),dist(o,o'))

[0060]

[0061] Where RD(o,o′) is the reachable distance from point o to point o′, core_dist(o′) is the core distance of point o′, and dist(o,o′) is the Euclidean distance between point o and point o′; σ B 2 (t) is the inter-class variance under threshold t, ω1(t) is the pixel ratio of the foreground class corresponding to threshold t, and ω2(t) is the pixel ratio of the background class corresponding to threshold t, μ1(t) is the mean of the foreground, μ2(t) is the mean of the background, μ T is the overall mean of the point cloud.

[0062] Finally, Gaussian Low-Pass Filtering (GLPF) is introduced to smooth the elevation sequence of the remaining photon points, further suppressing local random fluctuations and obtaining stable and reliable tidal flat prior elevation information. The calculation formula for the kernel density in the Gaussian low-pass filter is as follows:

[0063]

[0064] in, is the position h i The kernel density estimate at h i is the elevation value of the i-th photon point, h j is the elevation value of the j-th photon point, n is the total number of photons involved in the estimation, σ is the bandwidth of the kernel function (i.e., standard deviation), which controls the smoothness of the kernel density curve, and e is an exponential function.

[0065] Step 4: Construction and application of tidal flat terrain inversion model. The specific method is to first establish a grid with a size of 30m, and calculate the average value of ICESat-2 photon elevation within the grid as the prior elevation of the grid. Then, the prior elevation and its corresponding flooding frequency data are divided into a training set and a validation set (at a ratio of 8:2), and the tidal flat terrain inversion model is constructed using linear and nonlinear regression respectively. Figure 4 As shown, the constructed linear model is H l =-3.769f+2.957, R 2 is 0.787, and the nonlinear model is H nl =1.178f 3 +0.083f 2 –4.389f+3.065, R 2 Finally, the constructed terrain inversion model was used to convert the tidal flat inundation frequency in 2024 into tidal flat topography, and the tidal flat topography data in 2024 was obtained. Figure 5 As shown, Figure 5 (a) is the tidal flat topography result based on the linear model. Figure 5 (b) shows the tidal flat topography result based on the nonlinear model. Both results are in raster data format, with attribute values ​​representing tidal flat elevations. The elevation datum is EGM2008, with a spatial resolution of 30 meters and a temporal resolution of 1 year. The tidal flat topography inverted by this invention exhibits an overall spatial pattern of "low surroundings and high in the middle," which is consistent with the actual tidal flat topography and demonstrates the rationality and effectiveness of this invention.

[0066] Step 5: Accuracy evaluation. The accuracy of the tidal flat topography obtained by regression modeling is evaluated based on the airborne LiDAR data collected in 2024 and the ICESat-2 satellite altimeter data in the validation set. The calculation formula for accuracy evaluation is:

[0067]

[0068] Among them, RMSE is the root mean square error, MAE is the mean absolute error, Bias is the bias, R 2 is the coefficient of determination, y est Reverse the terrain elevation value for regression modeling, y obs is the elevation value of the airborne LiDAR data, is the average elevation value of the airborne LiDAR data, and n is the number of samples.

[0069] like Figure 6 As shown, Figure 6 (a) and Figure 6 (c) is the accuracy evaluation result of the linear model. Figure 6 (b) and Figure 6 (d) is the accuracy evaluation result of the nonlinear model. The results show that on the validation set, the RMSE of the linear model and the nonlinear model are 0.36m and 0.35m respectively, the MAE are 0.26m and 0.26m respectively, the Bias are 0.02m and 0.02m, and the R 2 Compared with the measured terrain, the RMSE of the linear model and the nonlinear model were 0.31m and 0.32m, the MAE were 0.21m and 0.22m, the Bias were 0.10 and 0.10m, and the R 2 The values ​​of 0.86 and 0.86 are respectively. Overall, the terrain inversion results based on the linear model are better than those based on the nonlinear model.

[0070] Example 2

[0071] like Figure 7 As shown in the figure, a tidal flat topography inversion system based on the collaboration of multi-source remote sensing and regression modeling includes:

[0072] The inundation frequency calculation module is used to pre-process multi-source remote sensing images, combine the land-water segmentation algorithm to map the tidal flats, and obtain the inundation frequency data of the tidal flats;

[0073] The tidal flat elevation extraction module is used to preprocess and denoise the lidar satellite altimeter data to obtain the prior elevation information of the tidal flat;

[0074] The tidal flat topography inversion module is used to construct linear and nonlinear tidal flat topography inversion models based on regression modeling, and to achieve accurate inversion of tidal flat topography by combining tidal flat inundation frequency.

[0075] Preprocessing of multi-source remote sensing images includes radiometric calibration, atmospheric correction, geometric precision correction, and unification of spatial resolution and coordinate system. Multi-source remote sensing images include Landsat-8 / 9, Sentinel-2, HJ-2A / B, and GF-1 optical images. Tidal flat inundation frequency is calculated using the following formula:

[0076]

[0077] Among them, F I is the flooding frequency of a pixel on the tidal flat, N tidal_flat is the frequency of a certain pixel being observed as a tidal flat, N sea_water is the frequency of observing seawater at a certain pixel on the tidal flat.

[0078] The lidar satellite altimeter data is from the ICESat-2 satellite altimeter. Preprocessing includes screening photon points with a confidence level of no less than 2 based on the confidence level, and spatial range screening based on the vector boundary of the target area. The elevation threshold for removing seawater photons is calculated according to the following formula:

[0079] T=H msl +2·σ msl

[0080] Among them, T is the elevation threshold, H msl is the average sea level height of photons obtained from the elevation histogram statistics, σ msl is the standard deviation of seawater photon elevation.

[0081] The tidal flat elevation extraction module analyzes the photon topographic characteristics of the ICESat-2 satellite altimeter data and uses a multi-level denoising and fusion method to accurately extract the prior elevation information of the tidal flat.

[0082] The density clustering method based on the adaptive reachable distance of the Otsu method is used to perform preliminary denoising on the tidal flat photons. The calculation formulas for the reachable distance and the Otsu method are as follows:

[0083] RD(o,o')=max(core_dist(o'),dist(o,o'))

[0084]

[0085] Where RD(o,o′) is the reachable distance from point o to point o′, core_dist(o′) is the core distance of point o′, and dist(o,o′) is the Euclidean distance between point o and point o′; σ B 2(t) is the inter-class variance under threshold t, ω1(t) is the pixel ratio of the foreground class corresponding to threshold t, and ω2(t) is the pixel ratio of the background class corresponding to threshold t, μ1(t) is the mean of the foreground, μ2(t) is the mean of the background, μ T is the overall mean of the point cloud.

[0086] The calculation formula of kernel density in Gaussian low-pass filtering of prior elevation information is as follows:

[0087]

[0088] in, is the position h i The kernel density estimate at h i is the elevation value of the i-th photon point, h j is the elevation value of the j-th photon point, n is the total number of photons involved in the estimation, σ is the standard deviation of the kernel function, which controls the smoothness of the kernel density curve, and e is an exponential function.

[0089] The general model of tidal flat topography inversion driven by regression modeling of the tidal flat topography inversion module is

[0090] H l =wf+b

[0091]

[0092] Among them, H l is a linear model, H nl is a nonlinear model, f is the flooding frequency, q is the highest order of the nonlinear model, w, b, u i , u0 is the model coefficient.

Claims

1. A tidal flat topography inversion method, characterized in that: The following steps are involved: Step 1: Acquire multi-source remote sensing images and ICESat-2 satellite altimeter data under light or no cloud conditions and perform preprocessing. Step 2: Use the support vector machine algorithm to classify the pre-processed multi-source remote sensing images and extract the flooding frequency information of the tidal flat; Step 3: By analyzing the photon topographic characteristics of ICESat-2 satellite altimeter data, a multi-level denoising and fusion method is used to accurately extract the prior elevation information of the tidal flat; Step 4: The tidal flat inundation frequency is mapped to the prior elevation information one by one, and a general tidal flat topography inversion model is constructed to map the inundation frequency to the tidal flat topography. Step 5: Use high-precision airborne LiDAR data to evaluate the accuracy of the inverted tidal flat topography.

2. A tidal flat topography inversion method according to claim 1, characterized in that: In step 1, the multi-source remote sensing images include Landsat-8 / 9, Sentinel-2, HJ-2A / B, and GF-1 optical images. The preprocessing of the multi-source remote sensing images includes radiometric calibration, atmospheric correction, geometric precision correction, and unification of spatial resolution and coordinate system. The preprocessing of ICESat-2 satellite altimeter data includes confidence screening and spatial range screening.

3. The tidal flat topography inversion method according to claim 1, characterized in that: In step 2, the tidal flat inundation frequency is calculated according to the following formula: F I =N sea_water / (N tidal_flat +N sea_water ) Among them, F I is the flooding frequency of a pixel on the tidal flat, N tidal_flat is the frequency of a certain pixel being observed as a tidal flat, N sea_water is the frequency of observing seawater at a certain pixel on the tidal flat.

4. The tidal flat topography inversion method according to claim 1, characterized in that: In step 3, the elevation threshold for removing seawater photons is calculated according to the following formula: T=H msl +2·s msl Among them, T is the elevation threshold, H msl is the average sea level height of photons obtained from the elevation histogram statistics, σ msl is the standard deviation of seawater photon elevation.

5. The tidal flat topography inversion method according to claim 1, characterized in that: In step 3, a density clustering method based on the adaptive reachable distance of the Otsu method is used to perform preliminary denoising on the tidal flat photons, wherein the calculation formulas of the reachable distance and the Otsu method are as follows: RD(o,o')=max(core_dist(o'),dist(o,o')) Where RD(o,o′) is the reachable distance from point o to point o′, core_dist(o′) is the core distance of point o′, and dist(o,o′) is the Euclidean distance between point o and point o′. is the inter-class variance under threshold t, ω1(t) is the pixel ratio of the foreground class corresponding to threshold t, and ω2(t) is the pixel ratio of the background class corresponding to threshold t, μ1(t) is the mean of the foreground, μ2(t) is the mean of the background, μ T is the overall mean of the point cloud.

6. The tidal flat topography inversion method according to claim 1, characterized in that: In step 3, the calculation formula of the kernel density in the Gaussian low-pass filter is as follows: in, is the position h i The kernel density estimate at h i is the elevation value of the i-th photon point, h j is the elevation value of the j-th photon point, n is the total number of photons involved in the estimation, and σ is the standard deviation of the kernel function, which controls the smoothness of the kernel density curve.

7. The tidal flat topography inversion method according to claim 1, characterized in that: In the step 4, a grid is established based on the pixel size of the tidal flat inundation frequency, and the prior elevation is averaged using grid filtering to obtain aligned prior elevation and inundation frequency data.

8. The tidal flat topography inversion method according to claim 1, characterized in that: In step 4, the general model for tidal flat topography inversion is as follows: H l =wf+b Among them, H l is a linear model, H nl is a nonlinear model, f is the flooding frequency, q is the highest order of the nonlinear model, w, b, u i , u0 is the model coefficient.

9. The tidal flat topography inversion method according to claim 1, characterized in that: In step 5, the calculation formula for accuracy evaluation is: Among them, RMSE is the root mean square error, MAE is the mean absolute error, Bias is the bias, R 2 is the coefficient of determination, y est Reverse the terrain elevation value for regression modeling, y obs is the elevation value of the airborne LiDAR data, is the average elevation value of the airborne LiDAR data, and n is the number of samples.

10. A tidal flat topography inversion system, characterized in that: include: The inundation frequency calculation module is used to pre-process multi-source remote sensing images, combine the land-water segmentation algorithm to map the tidal flats, and obtain the inundation frequency data of the tidal flats; The tidal flat elevation extraction module is used to preprocess and denoise the lidar satellite altimeter data to obtain the prior elevation information of the tidal flat; The tidal flat topography inversion module is used to construct linear and nonlinear tidal flat topography inversion models based on regression modeling, and to achieve accurate inversion of tidal flat topography by combining tidal flat inundation frequency.

Images