An improved gravity geological method for restoring seafloor topography at the sea-land boundary

Abstract

The present invention relates to the technical field of seabed topography inversion, and discloses an improved gravity geological method for restoring seabed topography in the sea-land boundary area, which utilizes the approximate linear relationship between short-wave gravity anomaly and sea depth, takes the external verification difference statistical result as a reference, and continuously adjusts the density difference between seawater and crust in the entire region, so that the sea depth and short-wave gravity anomaly reach a balance, thereby obtaining the optimal density difference constant of the entire irregular-shaped sea area; then, based on the functional relationship between short-wave gravity anomaly and sea depth, the seabed topography of the irregular-shaped sea area at the sea-land boundary area is inverted. The present invention is of great benefit to giving full play to the advantages of gravity data in acquiring seabed topography data in the sea-land boundary area, and accurately constructing a global seabed topography model.

Description

Technical Field

[0001] The invention relates to the technical field of seabed topography inversion, and in particular to an improved gravity geological method for restoring seabed topography in a sea-land boundary area. Background Art

[0002] At present, gravity inversion of seafloor topography has become the main means of constructing global seafloor topography models. Most of the global seafloor topography data released at home and abroad rely on gravity data recovery, such as STO_IEU2020, BAT_WHU2020, Sandwell series models, GEOCO series models, etc. The research on gravity inversion of seafloor topography mainly focuses on the restoration of seafloor topography in regular rectangular sea areas with full ocean coverage, and rarely involves the restoration of irregular seafloor topography at the junction of land and sea. The reasons may be that: 1) the admittance function method and the SA (Simulated Annealing) method generally process gravity data and topographic data in the frequency domain, involving Fourier forward and inverse transform operations, and are therefore usually applied to the inversion of seafloor topography in regular rectangular sea areas, but the restoration effect of irregular seafloor topography in the sea-land boundary area is not good; 2) the implementation of the seafloor topography inversion method generally involves the data gridding process, and the gridding results of the marine area data in the sea-land boundary area are bound to be interfered by the land area data, resulting in the weak robustness of the conventional seafloor topography inversion method to restore the marine part of the land-sea boundary area. However, on the one hand, the seafloor topography in the sea-land boundary area is an important part of the global seafloor topography, and is an important content that cannot be avoided in the independent establishment of the global seafloor topography model; on the other hand, with the continuous accumulation and enrichment of multi-source gravity data in the adjacent sea area and the continuous advancement of data processing technology, the quality of the available data related to the marine part of the sea-land boundary area will be greatly improved. It is foreseeable that the development and utilization of marine topography related data in the coastal area is imperative. Summary of the invention

[0003] In view of the shortcomings of gravity inversion algorithms for seabed topography in irregularly shaped sea areas at the sea-land boundary, such as difficulty in implementation and weak robustness, this invention proposes an improved gravity geological method for restoring the seabed topography in the sea-land boundary area based on the approximate linear relationship between seabed topography data and gravity data. This method is of great benefit for leveraging the advantages of gravity data in acquiring seabed topography data in the sea-land boundary area and accurately constructing a global seabed topography model.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] An improved gravity geological method for restoring seabed topography in the sea-land interface area, comprising:

[0006] Step 1: Divide the acoustic water depth measurement data into two categories according to the proportion. One category is the water depth discrete control point d k, used to participate in the construction of the seabed terrain model; the other type is the water depth discrete check point d j , used to evaluate the quality of seafloor topography models to obtain optimal density parameters;

[0007] Step 2: Set the initial value of the density difference constant △ρ0 based on d k , △ρ0 and reference water depth to solve the short-wave gravity anomaly at discrete control points

[0008] Step 3: Based on sea surface gravity anomaly data and Calculate long-wave gravity anomalies at discrete control points

[0009] Step 4: Gridding to obtain gridded long-wave gravity anomalies in the target area

[0010] Step 5: Based on Obtaining long wave gravity anomalies of the ocean grid in the land-sea boundary area

[0011] Step 6: Based on the sea surface gravity anomaly data and Obtain shortwave gravity anomalies of the ocean grid at the land-sea boundary

[0012] Step 7: Exploitation The inversion results of the gridded seafloor topography in the ocean part of the land-sea boundary area are obtained. gwo , that is, the seabed topography model of the land-sea interface;

[0013] Step 8: Based on d j , the accuracy of the seabed topography model in the land-sea interface area obtained by using root mean square evaluation;

[0014] Step 9: Adjust the density difference based on the set initial value of the density difference constant, and repeat steps 2 to 9 until the gridded seabed terrain model of the ocean part of the sea-land boundary area has the highest accuracy;

[0015] Step 10: Taking the density difference parameter corresponding to the minimum RMS value of the seabed topography model as the optimal density difference constant, the gridded seabed topography result of the corresponding sea-land boundary area is obtained in the same way as step 7;

[0016] Step 11: Use the nearest neighbor interpolation method to replace the numbers with values ​​greater than zero in the seabed topography results obtained in step 10 to obtain the final gridded seabed topography results of the marine part of the sea-land boundary area.

[0017] Furthermore, the step 2 comprises:

[0018] The shortwave gravity anomaly at the discrete control point is obtained using formula (1):

[0019]

[0020] In the formula, is the short-wave gravity anomaly of the control point; G is the gravitational constant; Δρ is the difference between the crust density and the seawater density at the study point; d k is the water depth at the control point; d is the reference water depth.

[0021] Furthermore, the step 3 comprises:

[0022] Step 3.1: Interpolate the sea surface gravity anomaly △g to the discrete control points to obtain the gravity anomaly △g at the control points k , and its interpolation formula is

[0023]

[0024] In the formula

[0025]

[0026] Where (x0, y0) represents the position of the discrete control point; x1, x2, y1 and y2 represent the position of the gravity anomaly model value point;

[0027] Step 3.2: Exploitation and △g k , the long wave gravity anomaly at the discrete control point is calculated according to the following formula:

[0028]

[0029] in, It is the control point long wave gravity anomaly.

[0030] Further, the step 6 comprises:

[0031] Step 6.1: Use the interpolation method described in step 3.1 to interpolate the sea surface gravity anomaly data in the land-sea boundary area to obtain the grid gravity anomaly △g of the ocean part of the land-sea boundary area. gwo ;

[0032] Step 6.2: Set △g gwo and The short-wave gravity anomaly of the ocean grid in the land-sea boundary area is obtained by subtraction.

[0033] Furthermore, the step 7 comprises:

[0034] The gridded seafloor topography result d of the ocean part of the sea-land boundary area is obtained by inversion according to the following formula gwo

[0035]

[0036] Where G is the gravitational constant and d is the reference water depth.

[0037] Furthermore, in step 8, the accuracy of the seabed terrain model at the land-sea boundary area is evaluated according to the following formula:

[0038]

[0039] Where d′ i represents the interpolation of the seabed terrain model at the i-th check point; Represents the sea depth value of the i-th check point.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] The present invention provides a strategy for modeling the seabed topography of irregular-shaped sea areas at the sea-land boundary. By using the approximate linear relationship between short-wave gravity anomaly and sea depth, and taking the external verification difference statistical results as a reference, the density difference between the seawater and the crust of the entire region is continuously adjusted, so that the sea depth and the short-wave gravity anomaly reach a balance, and then the optimal density difference constant of the entire irregular-shaped sea area is obtained; then, based on the functional relationship between the short-wave gravity anomaly and the sea depth, the seabed topography of irregular-shaped sea areas at the sea-land boundary is inverted.

[0042] The present invention has strong robustness, and the final inversion of the seabed topography of irregular-shaped sea areas at the sea-land boundary area is about 78% and 71% more accurate than that of SCS_DTU18 and SCS_ETOPO1, respectively. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A basic flow chart of an improved gravity geological method for restoring seabed topography at the sea-land boundary area according to an embodiment of the present invention;

[0044] Figure 2 This is a distribution diagram of water depth measurement points of a ship according to an embodiment of the present invention;

[0045] Figure 3 It is a distribution diagram of short-wave gravity anomaly and long-wave gravity anomaly at discrete control points in an embodiment of the present invention;

[0046] Figure 4 This is a long-wave gravity anomaly distribution map of the ocean part grid in the land-sea boundary area according to an embodiment of the present invention;

[0047] Figure 5 The distribution map of gravity anomaly and short-wave gravity anomaly in the ocean part of the land-sea boundary area according to the embodiment of the present invention;

[0048] Figure 6The final gridded seabed topography result SCS_iGGM of the ocean part of the sea-land boundary area obtained in the embodiment of the present invention;

[0049] Figure 7 These are example images of the DTU18 seafloor terrain model and the ETOPO1 seafloor terrain model. DETAILED DESCRIPTION

[0050] The present invention will be further explained below with reference to the accompanying drawings and specific embodiments:

[0051] like Figure 1 As shown, an improved gravity geological method for restoring the seabed topography in the sea-land interface area comprises:

[0052] Step 1: Divide the acoustic water depth measurement data into two categories according to a certain ratio; one is the water depth discrete control point d k , used to participate in the construction of the seabed terrain model; the other type is the water depth discrete check point d j , which is used to evaluate the quality of the seafloor topography model to obtain the optimal density parameters.

[0053] Step 2: Set the initial value of the density difference constant △ρ0 and solve the short-wave gravity anomaly at the discrete control point

[0054] The shortwave gravity anomaly at the discrete control point is obtained using formula (1):

[0055]

[0056] In the formula, is the short-wave gravity anomaly of the control point; G is the gravitational constant; Δρ is the difference between the crust density and the seawater density at the study point; d k is the water depth at the control point; d is the reference water depth, and the maximum water depth in the study area is usually selected as the reference water depth.

[0057] Step 3: Calculate the long-wave gravity anomaly at discrete control points

[0058] The step three comprises:

[0059] Step 3.1: Interpolate the sea surface gravity anomaly △g to the discrete control points to obtain the gravity anomaly △g at the control points k , and its interpolation formula is

[0060]

[0061] In the formula

[0062]

[0063] Where (x0, y0) represents the position of the discrete control point; x1, x2, y1 and y2 represent the position of the gravity anomaly model value point.

[0064] Step 3.2: Use step 2 to obtain the shortwave gravity anomaly △g of the control point s k and the gravity anomaly △g at the control point obtained in step 3.1 k , calculate the long wave gravity anomaly at the control point; the formula is as follows

[0065]

[0066] Among them, △g k It is the gravity anomaly of the sea surface control point; Long wave gravity anomaly for control points; It is the control point shortwave gravity anomaly.

[0067] Step 4: The long wave gravity anomaly at the control point in step 3 Gridding to obtain gridded long-wave gravity anomalies in the target area Specifically, the gridding method used is the "Surface" command of the GMT (Generic Mapping Tools) software to map the long-wave gravity anomaly at the control point in step 3. Gridding for target area gridding long wave gravity anomaly

[0068] Step 5: Use the interpolation method described in step 3.1 to grid the longwave gravity anomaly of the target area Interpolation to obtain long wave gravity anomalies of the ocean grid in the land-sea boundary area

[0069] Step 6: Obtain shortwave gravity anomalies of the ocean grid in the land-sea boundary area

[0070] The step six comprises:

[0071] Step 6.1: Use the interpolation method described in step 3.1 to interpolate the sea surface gravity anomaly △g of the target sea area to obtain the grid gravity anomaly △g of the sea part of the sea-land boundary area gwo .

[0072] Step 6.2: Substitute the △g obtained in step 6.1 gwo The obtained value in step 5 The short-wave gravity anomaly of the ocean grid in the land-sea boundary area is obtained by subtraction.

[0073] Step 7: Use the data obtained in step 6 The inversion results in the gridded seafloor topography of the ocean part of the land-sea boundary area. gwo; The calculation formula is as follows

[0074]

[0075] Wherein, the parameters have the same meaning as above.

[0076] Step 8: Evaluate the accuracy of the seabed terrain model at the land-sea boundary constructed in step 7. The formula is as follows

[0077]

[0078] Where d′ i represents the interpolation of the seabed terrain model at the i-th check point; Represents the sea depth value of the i-th check point.

[0079] Step 9: Adjust the density difference with a certain step size based on the set initial value of the density difference constant, and repeat steps 2 to 9 until the gridded seabed terrain model of the marine part of the sea-land boundary area is constructed with the highest accuracy.

[0080] Step 10: Take the density difference parameter corresponding to the minimum RMS value of the seabed topography model as the optimal density difference constant, and use equation (5) to invert and obtain the gridded seabed topography results of the marine part of the sea-land boundary area.

[0081] Step 11: Use the nearest neighbor interpolation / estimation method to replace the numbers with values ​​greater than zero in the seabed topography results obtained in step 10, and obtain the final gridded seabed topography results of the marine part of the sea-land boundary area. Specifically, the nearest neighbor interpolation / estimation method uses the coordinates of the known points and the coordinates of the points to be found to calculate the spatial distance between the known points and the points to be found, selects the sea depths of the 5 known points closest to the points to be found as input data, and interpolates / estimates the sea depth value of the points to be found by the distance weighted method.

[0082] As a specific implementable method, a 2°×2° (6°N~8°N, 115°E~117°E) area at the land-sea boundary of the South China Sea was selected as the research object. The method of the present invention was used to integrate shipborne bathymetric data and grav_31.1 gravity anomaly data to invert and restore the seabed topography data of the land-sea boundary area.

[0083] Step 1: The acoustic water depth measurement data is divided into two categories according to a certain ratio. One category is water depth discrete control points, and the other category is water depth discrete check points. In the embodiment, the check points and control points are 2467 and 9871 respectively. The distribution of control points and check points is as follows: Figure 2 As shown, the grey dots are control points and the black triangles are check points.

[0084] Step 2: Using the control point data in step 1, the reference water depth of the test area is 2971m. Set the initial value of the density difference constant △ρ0 to 1.00g / cm3 , obtain the shortwave gravity anomaly at discrete control points and long-wave gravity anomalies Respectively Figure 3 (a) and Figure 3 (b) as shown.

[0085] Step 3: Grid the long-wave gravity anomaly of the control point in step 2 to obtain the gridded long-wave gravity anomaly of the target area. Then the target area is gridded into long-wave gravity anomalies Interpolation to obtain long wave gravity anomalies of the ocean grid in the land-sea boundary area The long wave gravity anomaly of the ocean grid in the land-sea boundary area is as follows Figure 4 shown.

[0086] Step 4: Compare the target sea area grav_31.1 gravity anomaly data (called SCS_grav_31.1) with the data from step 3. By subtraction, we can obtain the short-wave gravity anomaly of the ocean grid in the sea-land boundary area. Named SCS_short_grav. SCS_grav_31.1 and SCS_short_grav are as follows Figure 5 (a) and Figure 5 (b) as shown.

[0087] Step 5: Use the SCS_short_grav obtained in step 4 to invert the gridded seabed topography of the ocean part of the land-sea boundary area. The RMS value of the corresponding seabed topography model verification is 827.43m.

[0088] Step 6: 0.01 g / cm 3 Repeat steps 2 to 5 for step size adjustment of density difference until the gridded seafloor terrain model of the sea-land boundary area is constructed with the highest accuracy (corresponding to the smallest check RMS). Finally, the optimal density difference constant (density difference parameter corresponding to the smallest check RMS of the seafloor terrain model) is 2.01 g / cm 3 .

[0089] Step seven, using the optimal density difference constant obtained in step six, invert and obtain the seabed topography model of the sea area of ​​the embodiment.

[0090] Step 8: Use the nearest neighbor interpolation / estimation method to replace the numbers greater than zero in the seabed topography results obtained in step 7, and obtain the final gridded seabed topography results of the marine part of the sea-land boundary area, named SCS_iGGM. Figure 6 shown.

[0091] The widely used DTU18 seafloor topography model and ETOPO1 seafloor topography model are introduced to participate in the horizontal comparison of models. Figure 7 As shown in the figure, the three seafloor topography models (sea depth models) are all well matched, and the SCS_DTU18 model has four obvious topographic bulges in the northwest compared with other models.

[0092] The statistical results of the differences between SCS_iGGM, SCS_DTU18 and SCS_ETOPO1 are shown in Table 1.

[0093] Table 1 Difference statistics between different seafloor topography models (unit: m)

[0094]

[0095] According to Table 1, the correlation coefficients of SCS_iGGM, SCS_DTU18 and SCS_ETOPO1 all exceed 0.99, indicating that the seabed topography data of the sea-land boundary restored by the method of the present invention are in good agreement with the overall data. The performance of SCS_iGGM, SCS_DTU18 and SCS_ETOPO1 models was further compared and analyzed, and the external accuracy of the seabed topography model was evaluated by using 3086 external water depth points that were not involved in the modeling. The above-mentioned sea depth model was interpolated to the measured water depth points, and the interpolation results were subtracted from the ship-measured sea depth values ​​at the measured water depth points. The corresponding difference statistical results are shown in Table 2.

[0096] Table 2 Statistical results of external conformity accuracy of the seabed terrain model (unit: m)

[0097]

[0098] It can be seen from Table 2 that the RMS of SCS_iGGM, SCS_DTU18, SCS_ETOPO1 and the external measured water depth values ​​are 40.35m, 183.79m and 138.15m respectively. It can be seen that the accuracy of SCS_iGGM is improved by about 78% and 71% respectively compared with SCS_DTU18 and SCS_ETOPO1, which further verifies the effectiveness of the method of the present invention.

[0099] The above is only a preferred embodiment of the present invention. It should be pointed out that a person skilled in the art can make several improvements and modifications without departing from the principle of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. An improved gravity geological method for restoring the seabed topography at the sea-land boundary, characterized in that: include: Step 1: Divide the acoustic water depth measurement data into two categories according to the proportion. One category is the water depth discrete control point d k , used to participate in the construction of the seabed terrain model; the other type is the water depth discrete check point d j , used to evaluate the quality of seafloor topography models to obtain optimal density parameters; Step 2: Set the initial value of the density difference constant △ρ0 based on d k , △ρ0 and reference water depth to solve the short-wave gravity anomaly at discrete control points The step 2 comprises: The shortwave gravity anomaly at the discrete control point is obtained using formula (1): In the formula, is the short-wave gravity anomaly of the control point; G is the gravitational constant; Δρ is the difference between the crust density and the seawater density at the study point; d k is the water depth at the control point; d is the reference water depth; Step 3: Based on sea surface gravity anomaly data and Calculate long-wave gravity anomalies at discrete control points The step 3 comprises: Step 3.1: Interpolate the sea surface gravity anomaly △g to the discrete control points to obtain the gravity anomaly △g at the control points k , and its interpolation formula is In the formula Where (x0, y0) represents the position of the discrete control point; x1, x2, y1 and y2 represent the position of the gravity anomaly model value point; Step 3.2: Exploitation and △g k , the long wave gravity anomaly at the discrete control point is calculated according to the following formula: in, Long wave gravity anomaly for control points; Step 4: Gridding to obtain gridded long-wave gravity anomalies in the target area Step 5: Based on Obtaining long wave gravity anomalies of the ocean grid in the land-sea boundary area Step 6: Based on the sea surface gravity anomaly data and Obtain shortwave gravity anomalies of the ocean grid at the land-sea boundary The step 6 comprises: Step 6.1: Use the interpolation method described in step 3.1 to interpolate the sea surface gravity anomaly data in the land-sea boundary area to obtain the grid gravity anomaly △g of the ocean part of the land-sea boundary area. gwo ; Step 6.2: Set △g gwo and The short-wave gravity anomaly of the ocean grid in the land-sea boundary area is obtained by subtraction. Step 7: Exploitation The inversion results of the gridded seafloor topography in the ocean part of the land-sea boundary area are obtained. gwo , that is, the seabed topography model of the land-sea interface; The step 7 comprises: The gridded seafloor topography result d of the ocean part of the sea-land boundary area is obtained by inversion according to the following formula gwo Where, G is the gravitational constant; d is the reference water depth; Step 8: Based on d j , the accuracy of the seabed topography model in the land-sea interface area obtained by using root mean square evaluation; In step 8, the accuracy of the seabed terrain model at the land-sea boundary area is evaluated according to the following formula: Where, d i ' represents the interpolation of the seabed terrain model at the i-th check point; represents the sea depth value of the i-th check point; Step 9: Adjust the density difference based on the set initial value of the density difference constant, and repeat steps 2 to 9 until the gridded seabed terrain model of the ocean part of the sea-land boundary area has the highest accuracy; Step 10: Taking the density difference parameter corresponding to the minimum RMS value of the seabed topography model as the optimal density difference constant, the gridded seabed topography result of the corresponding sea-land boundary area is obtained in the same way as step 7; Step 11: Use the nearest neighbor interpolation method to replace the numbers with values ​​greater than zero in the seabed topography results obtained in step 10 to obtain the final gridded seabed topography results of the marine part of the sea-land boundary area.

Images