Satellite altimetry and depth gauge combined seabed topography subsidence monitoring method

Abstract

The application discloses a kind of submarine topography subsidence monitoring methods of satellite altimetry and depth gauge combination, belong to marine surveying technical field, for submarine topography subsidence monitoring, including obtaining the time series satellite altimetry data and submarine depth gauge sounding data of monitoring sea area are preprocessed;Utilize tidal harmonic analysis method to calculate the amplitude and lag angle of astronomical component, eliminate periodic tidal fluctuation, obtain long-term average sea surface;The instantaneous sea surface height and long-term average sea surface height corresponding to the depth gauge position are obtained by bilinear interpolation;By calculating tidal deviation, the sounding data is uniformly reduced to long-term average sea surface datum;Submarine subsidence spatiotemporal sequence model is constructed, and submarine subsidence rate graph and cumulative subsidence amount surface are generated.The application is spatiotemporal matching and fusion by time series satellite altimetry data and depth gauge single-point high-precision data, realizes the vertical datum uniformity and high-precision transmission from sea surface to submarine, reduces the cost of traditional ship measurement or dense arrangement in-situ observation station.

Description

Technical Field

[0001] This invention discloses a method for monitoring seabed topographic subsidence by combining satellite altimetry and depth measurement, belonging to the field of marine measurement technology. Background Technology

[0002] Submarine subsidence is a critical process affecting global marine security and sustainable development. It is driven not only by natural factors such as plate tectonics, sediment compaction, and fluid extraction, but also increasingly significantly by human engineering activities such as oil and gas extraction, underground water injection, and large-scale coastal engineering projects. This slow and continuous subsidence alters seabed topography, potentially leading to the failure of subsea infrastructure, and poses a hidden and long-term threat by changing water levels, causing coastline erosion, wetland loss, and impacting flood control in coastal cities. Traditional seabed subsidence monitoring relies heavily on underwater bathymetry systems, which suffer from low accuracy, inefficiency, and high cost, making it difficult to obtain large-scale, high-precision overall deformation fields. With the development of satellite geodesy and marine geodesy technologies, revolutionary technical means have been provided for large-scale, continuous, and high-precision monitoring of seabed vertical displacement. This invention proposes a seabed topographic subsidence monitoring method that combines satellite altimetry and depth gauge. It can accurately acquire the vertical deformation of the seabed under a unified vertical reference, revealing key information on deep geological processes such as underground resource extraction and tectonic activity. It can also provide indispensable scientific basis for assessing the safety of seabed infrastructure and providing early warning of flood and erosion risks in coastal areas. Summary of the Invention

[0003] The purpose of this invention is to provide a method for monitoring seabed topographic subsidence by combining satellite altimetry and depth measurement, so as to solve the problems of low accuracy, poor efficiency, high cost and difficulty in obtaining large-scale overall deformation fields in existing seabed subsidence monitoring methods.

[0004] A method for monitoring seabed topographic subsidence by combining satellite altimetry and depth measurement includes:

[0005] S1. Determine the monitoring sea area and monitoring period, obtain the raw satellite altimetry data of the monitoring sea area, deploy an in-situ observation network using depth gauges, obtain the depth measurement data of the monitoring period, preprocess the raw satellite altimetry data and depth measurement data to obtain the time series data of sea surface height and the time series depth measurement data.

[0006] S2. Use the tidal harmonic analysis method to decompose the time series data of sea surface height, solve the signal decomposition results by the least squares method to obtain the amplitude and lag angle of the astronomical tide, and remove all periodic tidal fluctuations from the actual sea surface height based on the amplitude and lag angle to obtain the long-term mean sea surface height.

[0007] S3. Obtain the average sea level and instantaneous sea level at each depth gauge point through bilinear interpolation, and reduce the depth gauge bathymetry data to the long-term mean sea level by calculating the vertical offset.

[0008] S4. Construct a regular grid array for seabed subsidence to generate a seabed subsidence rate map and a cumulative subsidence surface for the monitored sea area.

[0009] S1 includes, S1.1, determining the monitoring sea area and monitoring period, and obtaining raw satellite altimetry data for the monitoring sea area. , For regular grid cells, , The number of regular grid cells;

[0010] S1 includes S1.2, which involves deploying an in-situ observation network based on the topography of the monitored sea area and satellite coverage. This network consists of depth gauges, and the depth measurement data collected by these depth gauges during the monitoring period is... , The serial number of the depth gauge. , For depth measurement quantity;

[0011] S1 includes, S1.3, and... Preprocessing is performed, including removing data points with confidence levels below a preset threshold, sequentially correcting for dry tropospheric delay, wet tropospheric delay, ionospheric delay, and sea state deviation, and eliminating orbital errors to obtain sea surface height time-series data. ;right Preprocessing is performed, including removing outliers, converting pressure values ​​to depth values, and performing sound velocity profile correction to obtain time-series bathymetry data. .

[0012] S2 includes, S2.1, for each within the monitoring area Through tidal harmonic analysis Perform signal decomposition:

[0013] ;

[0014] In the formula, for The long-term mean sea level height, The number of astronomical tides. For the index of tidal divisions, , for The The amplitude of each tidal component, for The The late angle of the tide, For observation time, for The astronomical phase angle of the tidal phase, for Measurement error.

[0015] S2 includes, S2.2, and Suppose have One observation value, for index, , establish containing and Tide height model function for each astronomical tide. :

[0016] ;

[0017] ;

[0018] In the formula, For tidal model parameters, It is the transpose symbol;

[0019] for Set initial parameter estimates :

[0020] ;

[0021] residual for:

[0022] ;

[0023] Calculate the deviation between the actual parameters and the initial parameters. :

[0024] ;

[0025] exist Perform a first-order Taylor expansion at this point:

[0026] ;

[0027] = ;

[0028] ;

[0029] ;

[0030] In the formula, , , To substitute variables;

[0031] Constructing the Jacobian matrix :

[0032] ;

[0033] and The observation equation is:

[0034] .

[0035] S2 includes, S2.3, and... Perform least squares calculation to construct the objective function. and make Minimum:

[0036] ;

[0037] In the formula, It is an L2 norm;

[0038] right Take the derivative and set it to 0 to construct the normal equation system matrix:

[0039] ;

[0040] The least squares solution is:

[0041] ;

[0042] Will As the initial parameter estimates for the iteration, based on Perform parameter estimation iterations, the first... The parameter estimates for the next iteration are:

[0043] ;

[0044] In the formula, For the first The parameter estimates for the next iteration. For the first The deviation of the next iteration;

[0045] Set iteration threshold ,when Stop iteration when the parameter estimate is obtained from the final iteration. calculate , and .

[0046] S3 includes, S3.1, determining the depth gauge location. coordinates , obtain The four nearest regular grid cells , , and , The coordinates are , The coordinates are , The coordinates are , The coordinates are ;

[0047] according to , , , and Calculate the weighted average of the distances. Mean sea level of the location and instantaneous sea level height :

[0048] ;

[0049] ;

[0050] In the formula, Mean sea level This refers to the instantaneous sea level height. The weighting coefficients are calculated as follows:

[0051] ;

[0052] In the formula, for The weight, for The weight, for The weight, for The weight.

[0053] S3 includes S3.2, representing the instantaneous sea level sequence as a vector based on the average sea level height and instantaneous sea level height of all depth gauge points:

[0054] ;

[0055] In the formula, This is a sequence vector of instantaneous sea level heights. For the first Instantaneous sea level height at each observation point For time point number, , The number of time-series observations at depth gauge points;

[0056] Represent the long-term mean sea level height as a vector:

[0057] ;

[0058] In the formula, This is a long-term mean sea level height sequence vector. For the first Long-term mean sea level height at each observation time point;

[0059] Represented by vector :

[0060] ;

[0061] In the formula, This is a time-series bathymetry data vector. For the first Time-series bathymetry data at each observation time point;

[0062] Tide level deviation sequence vector for:

[0063] ;

[0064] Substituting bathymetry data to the long-term mean sea level:

[0065] ;

[0066] In the formula, These are the depth values ​​of each depth gauge relative to the long-term sea level.

[0067] S4 includes, S4.1, calculating the predicted depth value of the predicted point:

[0068] ;

[0069] In the formula, For the prediction point, To predict the spatial location of the point, For the observation time of the prediction point, for The depth prediction value, These are the weighting coefficients. For the first The weighting coefficients of each depth measurement point, and , For the first depth gauge The depth value at a given time relative to the long-term sea level. For the first The spatial location of the depth gauge.

[0070] S4 includes S4.2, which defines the constrained optimization problem as finding the optimal weights by minimizing the variance of the prediction error while satisfying unbiasedness; the constrained optimization problem is transformed into solving a system of linear equations using the Lagrange multiplier method:

[0071] ;

[0072] In the formula, For spatiotemporal variability function, Let them be Lagrange multipliers; For row index, ,set up For column indexes, , Let be the spatiotemporal variability function between two depth gauge locations. This is the spatiotemporal variability function between depth gauge locations and prediction points; and The calculation method is as follows:

[0073] ;

[0074] ;

[0075] In the formula, For the first The spatial location of each depth gauge point For the first The spatial location of each depth gauge point For the first Observation time at each depth gauge point For the first Observation time at each depth gauge point;

[0076] Based on the principle of spatiotemporal kriging interpolation, using right Spatial interpolation and time series reconstruction are performed to form a regular grid array of seabed depth. By performing time difference calculation on the regular grid array of seabed depth, a regular grid array of seabed subsidence is generated. The regular grid array of seabed depth is used to monitor the seabed depth distribution at different times in the monitoring area, including spatial grid dimension and time series dimension. The regular grid array of seabed subsidence includes spatial grid dimension and time series dimension, and the array element value is the cumulative subsidence at the corresponding spatiotemporal point.

[0077] S4 includes S4.3, which generates a seabed subsidence rate map of the monitored sea area through time-series analysis based on a regular grid array of seabed subsidence data; within the monitoring period, any target time period is selected according to the coverage of the monitoring data. Through calculation Elevation difference, generating The cumulative settlement surface.

[0078] Compared to existing technologies, this invention offers the following advantages: By combining time-series satellite altimetry data with high-precision single-point depth gauge data through spatiotemporal matching and fusion, this invention effectively overcomes the limitations of single technologies, achieving a unified and high-precision vertical reference from the sea surface to the seabed. By utilizing satellite data as the primary source, only a limited number of depth gauges need to be deployed for calibration, significantly reducing the high costs and manpower required for traditional ship-based surveys or dense deployment of in-situ observation stations. Furthermore, by combining satellite altimetry data with depth gauge depth measurement data, this invention effectively separates sea surface change signals from seabed deformation signals, thereby enabling precise and direct monitoring of seabed subsidence. Attached Figure Description

[0079] Figure 1 A flowchart for monitoring seabed topographic subsidence by combining satellite altimetry and depth measurement;

[0080] Figure 2 Example of contour lines for cumulative settlement;

[0081] Figure 3 Example diagram of seabed settling rate;

[0082] Figure 4 Example diagram of the cumulative settlement surface. Detailed Implementation

[0083] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention are described clearly and completely below. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0084] A method for monitoring seabed topographic subsidence by combining satellite altimetry and depth measurement includes:

[0085] S1. Determine the monitoring sea area and monitoring period, obtain the raw satellite altimetry data of the monitoring sea area, deploy an in-situ observation network using depth gauges, obtain the depth measurement data of the monitoring period, preprocess the raw satellite altimetry data and depth measurement data to obtain the time series data of sea surface height and the time series depth measurement data.

[0086] S2. Use the tidal harmonic analysis method to decompose the time series data of sea surface height, solve the signal decomposition results by the least squares method to obtain the amplitude and lag angle of the astronomical tide, and remove all periodic tidal fluctuations from the actual sea surface height based on the amplitude and lag angle to obtain the long-term mean sea surface height.

[0087] S3. Obtain the average sea level and instantaneous sea level at each depth gauge point through bilinear interpolation, and reduce the depth gauge bathymetry data to the long-term mean sea level by calculating the vertical offset.

[0088] S4. Construct a regular grid array for seabed subsidence to generate a seabed subsidence rate map and a cumulative subsidence surface for the monitored sea area.

[0089] S1 includes, S1.1, determining the monitoring sea area and monitoring period, and obtaining raw satellite altimetry data for the monitoring sea area. , For regular grid cells, , The number of regular grid cells;

[0090] S1 includes S1.2, which involves deploying an in-situ observation network based on the topography of the monitored sea area and satellite coverage. This network consists of depth gauges, and the depth measurement data collected by these depth gauges during the monitoring period is... , The serial number of the depth gauge. , For depth measurement quantity;

[0091] S1 includes, S1.3, and... Preprocessing is performed, including removing data points with confidence levels below a preset threshold, sequentially correcting for dry tropospheric delay, wet tropospheric delay, ionospheric delay, and sea state deviation, and eliminating orbital errors to obtain sea surface height time-series data. ;right Preprocessing is performed, including removing outliers, converting pressure values ​​to depth values, and performing sound velocity profile correction to obtain time-series bathymetry data. .

[0092] S2 includes, S2.1, for each within the monitoring area Through tidal harmonic analysis Perform signal decomposition:

[0093] ;

[0094] In the formula, for The long-term mean sea level height, The number of astronomical tides. For the index of tidal divisions, , for The The amplitude of each tidal component, for The The late angle of the tide, For observation time, for The astronomical phase angle of the tidal phase, for Measurement error.

[0095] S2 includes, S2.2, and Suppose have One observation value, for index, , establish containing and Tide height model function for each astronomical tide. :

[0096] ;

[0097] ;

[0098] In the formula, For tidal model parameters, It is the transpose symbol;

[0099] for Set initial parameter estimates :

[0100] ;

[0101] residual for:

[0102] ;

[0103] Calculate the deviation between the actual parameters and the initial parameters. :

[0104] ;

[0105] exist Perform a first-order Taylor expansion at this point:

[0106] ;

[0107] = ;

[0108] ;

[0109] ;

[0110] In the formula, , , To substitute variables;

[0111] Constructing the Jacobian matrix :

[0112] ;

[0113] and The observation equation is:

[0114] .

[0115] S2 includes, S2.3, and... Perform least squares calculation to construct the objective function. and make Minimum:

[0116] ;

[0117] In the formula, It is an L2 norm;

[0118] right Take the derivative and set it to 0 to construct the normal equation system matrix:

[0119] ;

[0120] The least squares solution is:

[0121] ;

[0122] Will As the initial parameter estimates for the iteration, based on Perform parameter estimation iterations, the first... The parameter estimates for the next iteration are:

[0123] ;

[0124] In the formula, For the first The parameter estimates for the next iteration. For the first The deviation of the next iteration;

[0125] Set iteration threshold ,when Stop iteration when the parameter estimate is obtained from the final iteration. calculate , and .

[0126] S3 includes, S3.1, determining the depth gauge location. coordinates , obtain The four nearest regular grid cells , , and , The coordinates are , The coordinates are , The coordinates are , The coordinates are ;

[0127] according to , , , and Calculate the weighted average of the distances. Mean sea level of the location and instantaneous sea level height :

[0128] ;

[0129] ;

[0130] In the formula, Mean sea level This refers to the instantaneous sea level height. The weighting coefficients are calculated as follows:

[0131] ;

[0132] In the formula, for The weight, for The weight, for The weight, for The weight.

[0133] S3 includes S3.2, representing the instantaneous sea level sequence as a vector based on the average sea level height and instantaneous sea level height of all depth gauge points:

[0134] ;

[0135] In the formula, This is a sequence vector of instantaneous sea level heights. For the first Instantaneous sea level height at each observation point For time point number, , The number of time-series observations at depth gauge points;

[0136] Represent the long-term mean sea level height as a vector:

[0137] ;

[0138] In the formula, This is a long-term mean sea level height sequence vector. For the first Long-term mean sea level height at each observation time point;

[0139] Represented by vector :

[0140] ;

[0141] In the formula, This is a time-series bathymetry data vector. For the first Time-series bathymetry data at each observation time point;

[0142] Tide level deviation sequence vector for:

[0143] ;

[0144] Substituting bathymetry data to the long-term mean sea level:

[0145] ;

[0146] In the formula, These are the depth values ​​of each depth gauge relative to the long-term sea level.

[0147] S4 includes, S4.1, calculating the predicted depth value of the predicted point:

[0148] ;

[0149] In the formula, For the prediction point, To predict the spatial location of the point, For the observation time of the prediction point, for The depth prediction value, These are the weighting coefficients. For the first The weighting coefficients of each depth measurement point, and , For the first depth gauge The depth value at a given time relative to the long-term sea level. For the first The spatial location of the depth gauge.

[0150] S4 includes S4.2, which defines the constrained optimization problem as finding the optimal weights by minimizing the variance of the prediction error while satisfying unbiasedness; the constrained optimization problem is transformed into solving a system of linear equations using the Lagrange multiplier method:

[0151] ;

[0152] In the formula, For spatiotemporal variability function, Let them be Lagrange multipliers; For row index, ,set up For column indexes, , Let be the spatiotemporal variability function between two depth gauge locations. This is the spatiotemporal variability function between depth gauge locations and prediction points; and The calculation method is as follows:

[0153] ;

[0154] ;

[0155] In the formula, For the first The spatial location of each depth gauge point For the first The spatial location of each depth gauge point For the first Observation time at each depth gauge point For the first Observation time at each depth gauge point;

[0156] Based on the principle of spatiotemporal kriging interpolation, using right Spatial interpolation and time series reconstruction are performed to form a regular grid array of seabed depth. By performing time difference calculation on the regular grid array of seabed depth, a regular grid array of seabed subsidence is generated. The regular grid array of seabed depth is used to monitor the seabed depth distribution at different times in the monitoring area, including spatial grid dimension and time series dimension. The regular grid array of seabed subsidence includes spatial grid dimension and time series dimension, and the array element value is the cumulative subsidence at the corresponding spatiotemporal point.

[0157] S4 includes S4.3, which generates a seabed subsidence rate map of the monitored sea area through time-series analysis based on a regular grid array of seabed subsidence data; within the monitoring period, any target time period is selected according to the coverage of the monitoring data. Through calculation Elevation difference, generating The cumulative settlement surface.

[0158] The following description, in conjunction with the accompanying drawings, further illustrates the process of this invention. Figure 1 As shown, firstly, time-series satellite altimetry data and synchronous depth gauge bathymetry data are acquired, followed by data quality control and basic verification. Tidal harmonic analysis is performed on the sea surface height time-series data to obtain the amplitude and lag angle of each tidal constituent, and the long-term mean sea level is obtained. Interpolation is used to obtain the instantaneous sea level and long-term mean sea level corresponding to the depth gauge points, and the depth gauge bathymetry data is normalized to the long-term mean sea level. A spatiotemporal sequence model of seabed subsidence is constructed, and a seabed subsidence rate map of the monitored sea area is generated through time-series analysis. The cumulative subsidence surface is generated by calculating the elevation difference. Astronomical tidal constituents include the major lunar semi-diurnal tide (M2), the major solar semi-diurnal tide (S2), and the lunar-solar composite diurnal tide (K1), etc.

[0159] Figure 2 This paper presents a contour map of cumulative seafloor subsidence over a 5-year monitoring period generated by the method described in this invention, and sets... The preset confidence threshold is 95%. The monitoring area covers a range of 0-100 km eastward and 0-100 km northward, and is discretized using a 100×100 regular grid. The grayscale variations in the figure represent the cumulative seabed subsidence values ​​at different locations, in mm. The cumulative subsidence data is expressed as plane-filled contour lines, accurately quantifying the spatial distribution and gradient changes of subsidence. The subsidence center is located approximately 35 km eastward and 65 km northward, where the subsidence rate is highest, gradually decreasing outwards. The white triangles represent eight depth gauge locations, providing information such as monitoring duration, maximum and average subsidence, facilitating engineering assessments and risk zoning. Figure 3 A contour map of the spatial distribution of seafloor subsidence rates, generated according to the method described in this invention, is presented. The color mapping uses a grayscale scale, with light gray to dark black corresponding to subsidence rates ranging from 0 mm / year to -12 mm / year, visually reflecting the continuous spatial gradient change in subsidence rates. Figure 4 As shown, the cumulative settlement surface map visually presents the spatial distribution of cumulative settlement in the form of a topographic surface. The surface is processed with smooth interpolation, and the grayscale color mapping strictly corresponds to the surface elevation, which three-dimensionally shows the range, depth and transition relationship between the settlement area and the surrounding area.

[0160] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for monitoring the subsidence of the seabed topography by combining satellite altimetry and bathymetry, characterized in that, The method comprises the following steps: S1, determining a monitoring sea area and a monitoring period, obtaining satellite original height data of the monitoring sea area, deploying an in-situ observation network by using a depth gauge, obtaining sounding data of the monitoring period, pre-processing the satellite original height data and the sounding data to obtain sea surface height time series data and time series sounding data; S2, performing signal decomposition on the sea surface height time series data by using a tidal harmonic analysis method, solving the signal decomposition result by using a least square method, obtaining amplitudes and lag angles of astronomical component tides, removing all periodic tidal fluctuations from the actual sea surface height based on the amplitudes and the lag angles, and obtaining long-term average sea surface height; S3, obtaining average sea surface height and instantaneous sea surface height of each depth gauge point by using a bilinear interpolation, and calculating vertical offset to reduce the sounding data of the depth gauge to the long-term average sea surface height; S4, constructing a regular grid array of seabed subsidence, generating a seabed subsidence rate map and a cumulative subsidence amount surface of the monitoring sea area.

2. The method of claim 1, wherein the method is characterized by, S1 comprises S1.1, determining a monitoring sea area and a monitoring period, acquiring satellite raw altimetry data of the monitoring sea area , for a regular grid cell, , for a number of regular grid cells; S1 comprises S1.2, according to the monitoring sea area topography and satellite coverage, laying the in-situ observation network, the in-situ observation network is composed of depth gauge, the depth data of depth gauge in monitoring period is , is the serial number of depth gauge, , is the number of depth gauges; S1 includes S1.3, to The pretreatment includes eliminating data points with confidence lower than a preset threshold, sequentially performing dry troposphere delay correction, wet troposphere delay correction, ionosphere delay correction and sea state bias correction, and eliminating track errors to obtain sea surface height time series data ; to The pretreatment includes eliminating abnormal outliers, converting pressure values to depth values, and performing sound velocity profile correction to obtain time series depth data .

3. The method of claim 2, wherein the method is characterized by, S2 comprises, S2.1, for each of the monitoring regions signal decomposition by tidal harmonic analysis: signal decomposition by tidal harmonic analysis: ； In the formula, for The long-term mean sea level height The number of astronomical tides. For the index of tidal divisions, , for The The amplitude of each tidal component, for The The late angle of the tide, For observation time, for The astronomical phase angle of the tidal phase, for Measurement error.

4. The method of claim 3, wherein the method is characterized by, S2 comprises, S2.2, providing There are observations, are indices, , establishing a model function of tidal heights comprising and astronomical constituent : ； ； wherein is a tidal model parameter, is a transpose symbol; For Setting initial parameter estimates : ； residual of is: ； calculating a deviation of the real parameters from the initial parameters : ； At a first order Taylor expansion is made: ； = ； ； ； wherein , , are substitution variables; constructing the jacobian matrix : ； The observation equation for the observation of the position of the vehicle 1 is: The observation equation for the observation of the position of the vehicle 1 。 5. The method of claim 4, wherein the method further comprises, S2 comprises, S2.3, to perform a least square solution, construct an objective function and make min: ； wherein is the L2 norm; For Taking derivative and letting the derivative be 0, a normal equation matrix is constructed: ； The least squares solution is: ； Will As the initial parameter estimates for the iteration, based on Perform parameter estimation iterations, the first... The parameter estimates for the next iteration are: ； wherein is the parameter estimate for the th iteration, is the bias for the th iteration; Setting iteration threshold Stopping iteration when Stopping iteration when Calculating , and .

6. The method of claim 5, wherein the method is characterized by, S3 comprises, S3.1, determining the coordinates of the depth gauge point , , , , , , , , , , , , , , ; According to , , , and a weighted average of the distances is calculated the average sea surface height and the instantaneous sea surface height : ； ； wherein is the mean sea level height, is the instantaneous sea level height; is a weighting factor, calculated as: ； wherein is the weight of is the weight of is the weight of is the weight of is the weight of is the weight of is the weight of is the weight of 7. The method of claim 6, wherein the method further comprises, S3 comprises S3.2, representing the instantaneous sea surface height sequence by a vector according to the average sea surface height and the instantaneous sea surface height of all depth gauge points: ； wherein is the instantaneous sea surface height vector, is the instantaneous sea surface height at the th observation time point, is the time point index, , is the number of time series observations at the depth gauge site. representing the long-term average sea surface height by a vector: ； wherein is a long-term mean sea level height sequence vector, is the long-term mean sea level height at the observation time point; In vector notation : ； In the formula, is a time series sounding data vector, is the time series sounding data of the th observation time point. Tidal level bias sequence vector is: ； reducing the sounding data to the long-term average sea surface height: ； wherein are the depth values of each depth gauge relative to the long-term sea level.

8. The method of claim 7, wherein the method further comprises, S4 comprises S4.1, calculating a depth prediction value of a prediction point: ； wherein, is a prediction point, is a spatial position of the prediction point, is an observation time of the prediction point, is a depth prediction value of is a weight coefficient, is a weight coefficient of the th depth gauge point, and , is a weight coefficient of the th depth gauge is a depth value of the th depth gauge at the time instant relative to the long-term sea level, is a spatial position of the th depth gauge.

9. The method of claim 8, wherein the method further comprises, S4 comprises S4.2, defining a constraint optimization problem as solving optimal weights by minimizing the variance of prediction errors under the condition of satisfying unbiasedness; and converting the constraint optimization problem into solving a linear equation set by using a Lagrange multiplier method. ； wherein is the spatio-temporal variogram, is the Lagrange multiplier; let is the row index, , let is the column index, , is the spatio-temporal variogram between two depth gauge points, is the spatio-temporal variogram between a depth gauge point and a prediction point; and are calculated as follows: ； ； wherein is the spatial position of the th depth gauge point, is the spatial position of the th depth gauge point, is the observation time of the th depth gauge point, is the observation time of the th depth gauge point; Based on the principle of spatio-temporal Kriging interpolation, using spatial interpolation and time series reconstruction on a regular grid array of seafloor depth, through time difference calculation on the regular grid array of seafloor depth, a regular grid array of seafloor subsidence is generated; the regular grid array of seafloor depth is the distribution of seafloor depth at different times in the monitoring sea area, including spatial grid dimension and time series dimension; the regular grid array of seafloor subsidence contains spatial grid dimension and time series dimension, and the array element value is the cumulative subsidence of the corresponding spatio-temporal point.

10. The method of claim 9, wherein the method further comprises, S4 comprises S4.3, based on the subsea settlement gauge rule grid array, the subsea settlement rate map of the monitoring sea area is generated by time series analysis; in the monitoring period, according to the monitoring data coverage range, any target time period is selected , by calculating The elevation difference is generated The cumulative settlement surface.

Images