A method for restoring the thickness of passive continental margin initial crust

Abstract

The present application relates to the technical fields of marine geophysics, passive continental margin basin analysis and geodynamics numerical simulation, and discloses a passive continental margin initial crust thickness recovery method, comprising the following steps: obtaining seismic reflection profile interpretation, drilling data, heat flow data and tectonic zone partitioning and other multi-source data; establishing a present-day crustal structure model and a stretching-thinning-thermal evolution-interface burial depth response coupled numerical model; using decompaction correction and tectonic interface geometric parameters to establish initial crust thickness constraints; setting parameter ranges for different tectonic zone partitions; performing forward modeling with the initial crust thickness as the to-be-inverted parameter to obtain present-day simulation results; constructing an objective function according to the error between the simulation results and the measured constraints; updating the parameters using an optimization iterative algorithm until the convergence condition is met; and outputting the northern South China Sea continental margin initial crust thickness recovery results. The present application can reduce single-source constraint errors and improve the reliability of initial crust thickness recovery for strongly stretched continental margins.

Description

Technical Field

[0001] This invention relates to the fields of marine geophysics, passive continental margin basin analysis and geodynamic numerical simulation technology, and specifically to a method for restoring the initial crustal thickness of a passive continental margin. Background Technology

[0002] The northern continental margin of the South China Sea is a typical passive continental margin, having undergone multiple phases of crustal extension, thinning, and ocean basin expansion since the Cenozoic era. Its evolution not only constrains the spatiotemporal evolution of crustal thickness and Moho discontinuity depth but also controls the basin's sedimentary infilling process and thermal history. For marine oil and gas exploration, the northern continental margin of the South China Sea possesses a high geothermal flow background, profoundly influencing the hydrocarbon source rock maturation, migration, and accumulation processes in this region. Therefore, for oil and gas exploration in the northern continental margin of the South China Sea, constructing a reliable basin structure-thermal evolution model is crucial for identifying favorable source-reservoir-seal combinations, assessing the hydrocarbon generation-migration-accumulation sequence, and reducing drilling risks. Initial crustal thickness is a fundamental parameter for studying the mechanism of continental margin extension and thinning, basin formation and evolution, and lithospheric thermal evolution. The results of initial crustal thickness reconstruction directly affect the understanding of the extension coefficient, crustal thinning amount, and thermal history. Therefore, conducting initial crustal thickness reconstruction in the northern continental margin of the South China Sea has significant theoretical and practical value.

[0003] In existing technologies, initial crustal thickness recovery typically employs methods such as empirical assignment, seismic reflection profile interpretation, gravity inversion, and multi-parameter optimization inversion. Among these, optimization methods can establish geological or thermal evolution mathematical models and use measured parameters to invert unknown thicknesses. However, existing methods still have the following shortcomings: 1. Recovery results under single geophysical data constraints exhibit strong ambiguity; 2. Insufficient consideration of the heterogeneity and regional tectonic differences of the northern South China Sea margin; 3. Some recovery methods only focus on geometric thinning relationships, failing to fully couple thermal evolution, sedimentary load, and isostatic compensation processes; 4. A standardized, streamlined initial crustal thickness recovery technology suitable for the complex tectonic background of the northern South China Sea margin is lacking.

[0004] Therefore, it is necessary to propose a method for restoring the initial crustal thickness of the northern continental margin of the South China Sea based on multi-source data constraints, coupled numerical simulation and parameter inversion, so as to improve the stability, accuracy and engineering feasibility of the restoration results. Summary of the Invention

[0005] The purpose of this invention is to provide a passive method for restoring the initial crustal thickness of continental margins, which solves the problems of single constraints, multiple solutions, and insufficient regional adaptability in the restoration of initial crustal thickness in existing technologies. By introducing multi-source constraints such as seismic interpretation, heat flow, sedimentary thickness, and tectonic zone zoning, a coupled numerical model of extensional thinning-thermal evolution-interface burial depth response is established. An optimal iterative inversion strategy is adopted to achieve quantitative restoration of the initial crustal thickness of the northern South China Sea continental margin and improve the reliability and verifiability of the restoration of the initial crustal thickness of strongly extensional continental margins.

[0006] The technical problem solved by this invention is achieved by the following technical solution: A method for restoring the initial crustal thickness of a passive continental margin includes the following steps:

[0007] Step (1): Obtain multi-source geophysical and geological data, including seismic reflection profile interpretation data, present crustal thickness data, Moho discontinuity depth data, drilling data, magma data, heat flow data, sedimentary layer thickness data, basement depth data, and tectonic zone delineation data.

[0008] Step (2): Based on the multi-source geophysical and geological data, establish a current crustal structure model and a coupled numerical model of extensional thinning-thermal evolution-interface burial depth response for the study area;

[0009] Step (3): Based on the interpretation of seismic reflection profiles and drilling data, the depth of the same rift basement and the depth of the sedimentary bottom boundary after decompaction correction are used to establish the upper and lower limits of the initial crustal thickness.

[0010] Step (4): Divide the study area into tectonic zones and set the initial crustal structure and geophysical parameter ranges for different tectonic zones; Divide the study area into tectonic zones according to tectonic characteristics, from land to ocean: near-end zone, strong thinning zone, far-end zone, ocean-continent transition zone (COT), ocean basin; For the crustal structure and tectonic characteristics of different tectonic zones, set the initial crustal thickness, current crustal thickness, thermolithospheric thickness, current upper crustal thickness, current lower crustal thickness, and the range of values ​​for sedimentary layer, upper crust, and lower crustal density respectively; Among them, the near-end zone is represented by Enping Depression, the strong thinning zone is represented by Baiyun Depression and Baiyun East Depression is selected as a typical modeling area, the far-end zone is represented by Liwan Depression, and the study area also includes the outer edge highland area in the direction of ocean;

[0011] Step (5): Set the initial crustal thickness as the parameter to be inverted, and perform forward simulation based on the coupled numerical model to obtain the current simulation results; the model considers the influence of fault displacement, lithosphere deflection, sedimentary load and magma underplating on the crustal structure while performing forward simulation.

[0012] Step (6): Construct an objective function based on the difference between the current simulation results and the measured constraint data; the objective function is the weighted sum of at least two of the following: simulated crustal thickness error, simulated Moho depth error, simulated heat flow error, and simulated basement depth error.

[0013] Step (7): Update the parameters to be inverted using an optimization iterative algorithm, and repeat steps (5) to (6) until the preset convergence condition is met; the optimization iterative algorithm includes particle swarm optimization, gradient descent, quasi-Newton algorithm, genetic algorithm, simulated annealing algorithm or a combination thereof; the parameters to be inverted include at least one of the following: initial crustal thickness zoning parameters, fault extension parameters, or magma underplating thickness parameters.

[0014] Step (8): Output the initial crustal thickness recovery results of the northern continental margin of the South China Sea.

[0015] The specific steps in step (2) above for constructing the current crustal structure model are as follows: First, based on step (1), according to the tectonic trend of the study area, the distribution direction of the seismic reflection profile, and the research objectives, a two-dimensional research profile passing through the Baiyun East Depression, the Liwan Depression, and the ocean-facing outer high region of the study area is selected as the modeling profile; a two-dimensional rectangular coordinate system is established along the modeling profile, wherein the horizontal direction is... Direction, vertical depth direction is Direction, among which Corresponding to mean sea level, The area is the air layer. The area consists of the seawater layer and the underlying geological bodies;

[0016] Based on seismic reflection profile interpretation data, well layer data, current crustal thickness data, Moho depth data, and basement depth data, the main geological interface control points on the modeling profile are extracted. The main geological interfaces include at least the sea level, seafloor interface, sedimentary basement interface, Conrad discontinuity, Moho discontinuity, and lithospheric boundary. The current crustal structure model is divided into eight layers from top to bottom: air layer, seawater layer, sedimentary layer, upper crust, lower crust, lower crustal high-velocity body, lithospheric mantle, and asthenosphere.

[0017] The specific steps in step (2) above for constructing the coupled numerical model of stretching thinning—thermal evolution—interface burial depth response are as follows:

[0018] Based on the existing crustal structure model, a coupled numerical model of extensional thinning, thermal evolution, and interface burial depth response is further established. This coupled numerical model includes a crustal extensional thinning sub-model, a lithospheric thermal evolution sub-model, a sedimentary load sub-model, and a deflection and isostatic compensation sub-model. By coupling these sub-models, a coupled numerical model of extensional thinning, thermal evolution, sedimentary load, and deflection and isostatic compensation is formed. Specifically, the existing crustal structure model provides the geometric framework and material property distribution; the crustal extensional thinning sub-model provides the crustal thinning amount and extension factor; the lithospheric thermal evolution sub-model provides the temperature field and heat flow response; the sedimentary load sub-model provides the sedimentary loading amount; the deflection and isostatic compensation sub-model calculates the burial depth response of the basement interface, Moho discontinuity, and lithospheric interface based on load changes and density structure; and the magma underplating response expression is used to describe the influence of high-velocity anomalies at the bottom of the lower crust on density structure, thermal structure, and interface response.

[0019] The inputs to the coupled numerical model include current crustal structure geometry data, material property parameters, fault interpretation data, sedimentary layer thickness data, heat flow data, magmatic activity data, and Moho discontinuity depth data; the outputs of the coupled numerical model include simulated crustal thickness, simulated basement depth, simulated Moho discontinuity depth, and spatial distribution of high-velocity bodies in the lower crust.

[0020] Specifically, in step (3) above: let the current thickness of the sedimentary layer be... The thickness of the deposited layer after decompaction Represented as:

[0021] ;

[0022] in, and These represent the thickness of the sedimentary layer after decompaction and the current thickness of the sedimentary layer, respectively, in units of... ; and These are the average porosity and initial porosity of the sedimentary layer, respectively, and are dimensionless. Furthermore, the average porosity is obtained statistically from drilling data, while the initial porosity is obtained by inversion using an exponentially decaying compaction model. The compaction model is as follows:

[0023] ;

[0024] in, For burial depth Porosity at a given location, dimensionless; The unit is ; This is the compaction coefficient, in units of... e is the natural constant; generally, the compaction coefficient Preferred concentration: 0.3-0.5 ;

[0025] Based on the interpretation of seismic reflection profiles and the calculation of the predicted burial depth of key interfaces at each period using drilling data, the burial depth of the same rift basement is used as the lower limit constraint for the initial crustal thickness inversion, and the burial depth of the sedimentary bottom boundary after decompaction correction is used as the upper limit constraint for the initial crustal thickness inversion.

[0026] Specifically, in step (6) above:

[0027] Let the parameter vector to be inverted be The objective function is The objective function is dimensionless and is a weighted sum of the simulated crustal thickness error and the simulated basement depth error.

[0028] ;

[0029] in, and These are the simulated crustal thickness and the measured crustal thickness, respectively, with units of... ; and These are the simulated and measured foundation depths, respectively, with units of [units not specified]. ; and These are weighting coefficients, dimensionless.

[0030] Let the first The particle in the first The velocity vector and position vector of the line are respectively and Its optimal position is The global optimal position is Since it is dimensionless, the update formula for the particle swarm optimization algorithm is:

[0031] ;

[0032] ;

[0033] in, For inertial weights, For individual learning factors, As a group learning factor, and A random number between 0 and 1, dimensionless;

[0034] When the particle swarm optimization algorithm iterates into the region where the objective function is relatively optimal, the parameters are further updated using the local gradient correction formula:

[0035] ;

[0036] in, and These represent the local optimal parameters at the th... generation and first The value of the substitute, Step size, The gradient of the objective function with respect to the parameters to be inverted is dimensionless.

[0037] During the iteration process, a global search is first performed using the particle swarm optimization algorithm. When the objective function decreases by less than a preset decrease threshold for several consecutive iterations, local gradient correction is initiated to refine the current optimal solution. This process is repeated until the change in the objective function is less than a preset convergence threshold or the preset maximum number of iterations is reached. The convergence condition is defined as the relative rate of change of the objective function in two adjacent iterations being less than a preset threshold, specifically expressed as follows:

[0038] ;

[0039] in, and They represent the first Second and third The objective function value of the next iteration. The threshold value is dimensionless.

[0040] When the above conditions are continuously less than the preset decreasing threshold When the above conditions are less than the preset convergence threshold, local gradient correction is initiated; when the above conditions are less than the preset convergence threshold, local gradient correction is initiated. Or reach the preset maximum number of iterations Stop iterating when the time is right. , All are dimensionless parameters and satisfy the following conditions: > > >0.

[0041] The output results in step (8) above include the initial crustal thickness distribution and its confidence interval of the northern continental margin of the South China Sea, as well as the distribution of tectonic zone extension coefficients; wherein, the confidence interval of the initial crustal thickness is used to characterize the uncertainty of the inversion results, and its calculation method is as follows: under the premise of satisfying the preset convergence conditions, retain the feasible solution set obtained by multiple inversions, statistically analyze the initial crustal thickness results of each region, and take its quantile interval as the confidence interval; while the distribution map of tectonic zone extension coefficients is used to characterize the crustal extension intensity of tectonic zones in different regions; specifically, the confidence interval is expressed as:

[0042] ;

[0043] in, For position The confidence interval for the initial crustal thickness. and These represent the lower and upper quantiles of the initial crustal thickness inversion result at this location, respectively. For the preset significance level, preferably, The value is typically 0.05.

[0044] Specifically, its elongation coefficient is expressed as:

[0045] ;

[0046] in, is the extension coefficient of the tectonic zone, dimensionless; This represents the initial crustal thickness. The current thickness of the Earth's crust, in units of 1 / 2. .

[0047] Compared with the prior art, this application has the following advantages:

[0048] 1. This application reduces the ambiguity caused by single data recovery by using multi-source geophysical and geological data as constraints.

[0049] 2. This application improves the geological rationality of the initial crustal thickness restoration results by coupling simulation of extensional thinning, thermal evolution, sedimentary load and isostatic compensation.

[0050] 3. This application uses a zonal parameter initialization method to adapt to the longitudinal variation characteristics of the northern continental margin of the South China Sea.

[0051] 4. This application improves the automation and repeatability of the recovery results through optimal iterative inversion.

[0052] 5. This application is applicable to the restoration of passive continental margin structures, the evaluation of marine oil and gas resource potential, and basin evolution and exploration studies. Attached Figure Description

[0053] Figure 1 This is a schematic flowchart of the passive continental margin initial crust thickness restoration method of the present invention;

[0054] Figure 2 A schematic diagram of the unified processing flow for multi-source data;

[0055] Figure 3 A schematic diagram illustrating a unified scale standard for multi-source data;

[0056] Figure 4 A unified gridded representation of the results from multi-source data in the study area;

[0057] Figure 5 This is a schematic diagram of the current Earth's crustal structure.

[0058] Figure 6 This is a simulation of the current crustal structure of the Baiyun East Depression-Liwan Depression in the Baiyun Depression based on the Underworld2 coupled numerical model.

[0059] Figure 7 This is a cross-sectional view of the current crustal structure of the Baiyun East Depression-Liwan Depression, based on the forward modeling of the Underworld2 coupled numerical model.

[0060] Figure 8 This is a cross-sectional view of the initial crustal structure of the Baiyun East Depression-Liwan Depression, which was inverted and restored based on the Underworld2 coupled numerical model.

[0061] Figure 9 This is a map showing the initial crustal thickness, confidence interval (95% CI), and extension coefficient β distribution of the Baiyun East Depression-Liwan Depression, which were inverted and restored based on the Underworld2 coupled numerical model. Detailed Implementation

[0062] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0063] Example 1: The flowchart of a passive continental margin initial crust thickness restoration method is as follows: Figure 1 As shown, the specific steps are as follows:

[0064] Step (1): Select typical two-dimensional seismic reflection interpretation profiles of the Enping Depression-Baiyun Depression-Liwan Depression and the outer edge high area in the ocean direction as the research object, and obtain multi-source geophysical and geological data of the area. The multi-source geophysical and geological data include seismic reflection profile interpretation data, current crustal thickness data, Moho discontinuity depth data, drilling data, magma data, heat flow data, sedimentary layer thickness data, basement burial depth data and tectonic zone division data.

[0065] Unified coordinate transformation, spatial registration, outlier removal, and scale standardization are performed on multi-source data to ensure consistency in spatial location and value range across different sources; the specific process is as follows: Figure 2 As shown.

[0066] Specifically, firstly, seismic reflection profile interpretation, structural zoning data, and related drilling data are uniformly converted to the same coordinate system, and spatial registration is performed according to the study area boundary, allowing various data to be overlaid and analyzed within the same spatial reference frame. Furthermore, outlier removal or correction is performed on data that significantly deviates from the normal distribution range, and scale standardization is applied to data with different value ranges. The processed results are as follows: Figure 3 As shown in Table 1, a unified scaling standard is used to eliminate the impact of the order-of-magnitude differences between parameters on subsequent inversion calculations.

[0067] Table 1. Attribute Comparison Table Before and After Multi-Source Data Processing

[0068]

[0069] Furthermore, the processed data is represented in a grid format, with grid cells serving as the unified computational unit; the grid spacing is set to 10 based on the study area scale and data resolution. -100 Any suitable value in, preferably 50 Through the above processing, a unified gridded representation of the multi-source data for the study area is formed, providing an input basis for subsequent forward and inverse retrievals. In this embodiment, the depth benchmark is the mean sea level, and the vertical coordinates are positive downwards. The results are as follows: Figure 4 As shown.

[0070] Step (2): Based on the multi-source geophysical and geological data, establish a current crustal structure model and a coupled numerical model of extensional thinning-thermal evolution-interface burial depth response for the study area;

[0071] First, based on step (1), according to the structural orientation of the study area, the distribution direction of the seismic reflection profile, and the research objectives, a two-dimensional research profile passing through the Baiyun East Depression, the Liwan Depression, and the high area on the outer edge of the study area in the ocean direction is selected as the modeling profile. A two-dimensional rectangular coordinate system is established along the modeling profile, wherein the horizontal direction is... Direction, vertical depth direction is Direction: The horizontal length of the 2D model is set to 360. The vertical calculation depth is set to 150. ,in Corresponding to mean sea level, The area is the air layer. The area consists of the seawater layer and the underlying geological bodies.

[0072] Based on seismic reflection profile interpretation data, well layer data, current crustal thickness data, Moho discontinuity depth data, and basement depth data, key geological interface control points are extracted from the modeling profile. These key geological interfaces include at least the sea level, seafloor interface, sedimentary basement interface, Conrad discontinuity, Moho discontinuity, and lithospheric boundary. For any geological interface, its discrete interpretation points along the profile direction are represented as follows:

[0073] ;

[0074] in For discrete points on the profile, Indicates the first A geological interface, This indicates the first one on this interface. Each control point is dimensionless. To determine the position of the control point along the profile direction, For its corresponding depth, the unit is . .

[0075] Interpolation is performed on multiple control points on the same geological interface to obtain the continuous function of that interface:

[0076] ;

[0077] in, Indicates the first The geological interface is in a horizontal position The depth at that location It is its discrete interface The continuous function expression, with units of The interpolation methods include linear interpolation, piecewise linear interpolation, or spline interpolation. Preferably, since the seismic interpretation interface exhibits piecewise variations and the number of control points for some interfaces is limited, a piecewise linear interpolation method is used to construct a continuous interface.

[0078] Based on the aforementioned interface functions, the current crustal structure model of the study area is divided into eight layers from top to bottom: air layer, oceanic layer, sedimentary layer, upper crust, lower crust, lower crustal high-velocity body, lithospheric mantle, and asthenosphere. The model is as follows: Figure 5As shown. Specifically, the layer above sea level is the air layer; the layer between sea level and the seafloor interface is the seawater layer; the layer between the seafloor interface and the sedimentary basement interface is the sedimentary layer; the layer between the sedimentary basement interface and the Conrad discontinuity is the upper crust; the layer between the Conrad discontinuity and the Moho discontinuity is the lower crust; the layer between the Moho discontinuity and the lithospheric boundary is the lithospheric mantle; and the layer below the lithospheric boundary is the asthenosphere. The lower crustal high-velocity body is selectively set based on the interpretation results of seismic reflection profiles or magmatic activity data in the study area. When strong lower crustal reflection zones or evidence of magma underplating exist in the study profile, a lower crustal high-velocity body is set near the bottom of the lower crust; when no such evidence exists, no lower crustal high-velocity body is set, or the thickness of the lower crustal high-velocity body is set to zero. The spatial extent of the lower crustal high-velocity body is defined by its top and bottom interfaces:

[0079] ;

[0080] in, This is the function of the top interface of the high-velocity body in the lower crust. This represents the interface function of the high-velocity lower crust. This leads to the geometric framework of the current crustal structure model for the study area. It is important to note that... Figure 5 The eight-layered structure of the current crustal structure model is only used to illustrate the layer division method of the model and does not limit the specific burial depth of each layer interface; the actual burial depth of each layer interface is determined by seismic reflection profile interpretation data, drilling data, Moho depth data and basement burial depth data.

[0081] Secondly, the two-dimensional current crustal structure geometric framework is discretized into multiple grid cells. Taking the Baiyun East-Liwan Depression profile as an example, the model computational domain is set to 360°. ×150 Discretization is performed using a regular rectangular grid, with a grid cell size of approximately 50×20. For the regular rectangular grid, the horizontal grid spacing is expressed as:

[0082] ;

[0083] The vertical grid spacing is:

[0084] ;

[0085] in, and These represent the grid spacing in the horizontal and vertical directions, respectively, in units of [units missing]. ; and These represent the horizontal length and vertical depth of the model, respectively, in units of... ; This represents the number of grid cells in the horizontal direction. This represents the number of grid cells in the vertical direction, expressed in units of cells; the coordinates of the center point of any grid cell are represented as follows:

[0086] ;

[0087] in, and These represent the grid numbers in the horizontal and vertical directions, respectively.

[0088] It should be noted that the simulation is performed using the Underworld2 numerical simulation platform. Within the numerical simulation framework, the geometric and physical property information of the mesh cells is mapped to particle or particle swarm attributes. That is, when building the model using the Underworld2 platform, particles or particle swarms are placed within the mesh cells. These particles carry material identification numbers, spatial locations, and physical property parameters. The number of particles uses the default particle swarm placement method of the Underworld2 platform, but can also be set to multiple particles per mesh cell depending on the required computational accuracy. Preferably, 5-30 particles are placed per mesh cell. Each particle includes at least spatial location, layer attribute, material identification number, density attribute, thermal conductivity attribute, specific heat capacity attribute, radioactive heat generation rate attribute, and reference viscosity attribute.

[0089] For any particle Let its coordinates be Based on the spatial position of the particle relative to various geological interfaces, its material layer is determined and a material value is assigned. The specific discrimination rules are as follows:

[0090] when When located above sea level, the particles are assigned the characteristics of air layer materials;

[0091] when When located between the sea level and the seabed interface, the particles are assigned the value of seawater layer material;

[0092] when When located between the seabed interface and the sedimentary substrate interface, the particles are assigned the characteristics of the sedimentary layer material;

[0093] when When located between the sedimentary basement interface and the Conrad surface, the particles are assigned the characteristics of upper crustal material;

[0094] when When the particle is located between the Conrad surface and the Moho surface and is not within the range of the lower crust high-velocity body, the particle is assigned the value of lower crust material;

[0095] when When located between the top and bottom interfaces of the lower crustal high-velocity body, the particles are assigned the characteristics of the lower crustal high-velocity body material.

[0096] when When located between the Moho discontinuity and the lithospheric boundary, the particles are assigned the characteristics of lithospheric mantle material;

[0097] when When located below the lithospheric boundary, the particles are assigned as asthenospheric materials.

[0098] Using the above method, the stratigraphic geometry information obtained from seismic interpretation is converted into grid distribution and particle material distribution in the numerical model. After the above gridding and material assignment, the current seismic profile data of the Baiyun East-Liwan Depression in the study area is generated in the Underworld2 numerical simulation platform as follows: Figure 6 The image shows a current two-dimensional numerical simulation model of the Earth's crust.

[0099] The density and thermal parameters of each layer of material were assigned values ​​based on existing regional research results, rock physical properties, heat flow data, and seismic interpretation results. These physical properties include at least density, thermal conductivity, specific heat capacity, radioactive heat generation rate, and reference viscosity. In this embodiment, the typical ranges of physical properties for each layer of material are shown in Table 2.

[0100] Table 2. Typical material property ranges for each layer in the current crustal structure model.

[0101]

[0102] It should be noted that the parameters shown in Table 2 are typical parameter ranges used for establishing the numerical model in this embodiment. In specific calculations, the parameters are selected within these ranges based on existing rock property test results, regional heat flow data, and existing geodynamic research findings for the specific study area. The grid cell parameters listed in Table 2 correspond to the spatial location, stratigraphic affiliation, and physical properties of representative particles in the model, illustrating the differences in physical properties between different material stratigraphic layers. During the actual calculation process, all grid cells and particles are automatically assigned corresponding material properties according to their spatial location and stratigraphic affiliation.

[0103] Based on the existing crustal structure model, a coupled numerical model of extensional thinning, thermal evolution, and interface burial depth response is further established. This coupled numerical model includes a crustal extensional thinning sub-model, a lithospheric thermal evolution sub-model, a sedimentary load sub-model, and a deflection and isostatic compensation sub-model.

[0104] The crustal extension thinning sub-model is used to describe the thickness changes of the crust during rift extension. For any horizontal position on the cross-section... The current thickness of the Earth's crust is expressed as:

[0105] ;

[0106] in, The current thickness of the Earth's crust, The current depth of the Moho. The depth of the sedimentary basement is given in units of 1 / 2. .

[0107] The initial crustal thickness is expressed as... ,unit The crustal extension factor is then expressed as:

[0108] ;

[0109] The amount of crustal thinning is expressed as:

[0110] ;

[0111] in, is the crustal extension coefficient, dimensionless; In this step, the upper and lower limits of the variables used as input variables and variables to be updated in subsequent forward and inverse simulations are not determined.

[0112] Fault displacement data were obtained based on fault interpretation results from seismic reflection profiles of the study area. In this embodiment, faults are not directly plotted as independent material horizons in the Underworld2 model map, but rather used as geometric constraint inputs for the crustal extension thinning sub-model. Specifically, fault displacements, fault combination patterns, and fault-controlled sedimentary thickness variations identified on the seismic profiles are used to constrain local crustal extension intensity, sedimentary basement morphology, and Moho response. For faults whose displacement and dip angle can be quantitatively extracted, their horizontal extension is expressed as:

[0113] ;

[0114] in, For the first The horizontal extension of the fault. For the first The displacement of the fault is expressed in units of 1 / 2000. ; For the first The dip angle of the fault, in units of The cumulative horizontal extension caused by faults within a structural unit is expressed as:

[0115] ;

[0116] in, This represents the cumulative horizontal extension caused by faults within this structural unit, expressed in units of... ; This represents the number of faults involved in extension within the tectonic unit, expressed in units of faults. For faults whose displacement is difficult to extract quantitatively, fault development density, fault-controlled sedimentary thickness variations, and fault combination patterns are used as semi-quantitative constraints to limit the extensional thinning intensity at the corresponding locations.

[0117] The lithospheric thermal evolution sub-model is used to describe the evolution of the temperature field and heat flow under the background of rift extension and lithospheric thinning. The temperature field in the model satisfies the heat conduction governing equation:

[0118] ;

[0119] in, Temperature, unit: ; Time, in seconds; Density, unit: ; Specific heat capacity, unit: ; Thermal conductivity, in units of ; The rate of radioactive heat generation is expressed in units of [missing information]. ; For divergence operators, This is the gradient operator. Different thermal property parameters are used for different material layers: air layer, seawater layer, sedimentary layer, upper crust, lower crust, lower crust high-velocity body, lithospheric mantle, and asthenosphere are respectively assigned corresponding density, thermal conductivity, specific heat capacity, and radioactive heat generation rate.

[0120] The top and bottom of the model are set as fixed temperature boundaries; the left and right sides of the model are set as adiabatic boundaries or zero-horizontal heat flux boundaries. Preferably, the temperature at the top of the model is determined based on the seabed temperature or surface temperature, with a range of 0-25℃; the temperature at the bottom of the model is determined based on the asthenosphere temperature, with a range of 1250-1350℃; the left and right boundaries are zero-horizontal heat flux boundaries. The seabed heat flux is calculated from the temperature gradient.

[0121] ;

[0122] in, For position The heat flux value at the location, in units of ; For depth, the unit is Heat flow data is used to constrain the rationality of thermal conductivity, lithosphere thickness, radioactive heat generation rate, and thermal anomaly distribution in the model.

[0123] The sediment loading sub-model is used to describe the impact of sediment loading on the burial depth of the crustal interface. For any horizontal position on the profile... The thickness of the sedimentary layer is expressed as:

[0124] ;

[0125] in, For the thickness of the sedimentary layer, The depth of the seabed interface. The depth of the sedimentary basement is given in units of 1 / 2. .

[0126] The variation of sediment porosity with depth is expressed as:

[0127] ;

[0128] in, and Depth Porosity at the given location and initial porosity are both dimensionless. This is the compaction coefficient, in units of... The density of the sedimentary layer is expressed as:

[0129] ;

[0130] in, The density of the sedimentary layer, For the density of sedimentary particles, The density of seawater is given in units of 1 / 2. Deposition load is expressed as:

[0131] ;

[0132] Among them, deposition load The unit is ; The depth of the seabed interface. The depth of the sedimentary basement is given in units of 1 / 2. ; This refers to the density of the sedimentary layer, in units of... ; Let be the acceleration due to gravity, taken as 9.8. The depositional load, as a load term in the interface burial depth response calculation, is used to calculate the basement interface and Moho surface response caused by deposition. It is important to note that during the calculation, the interface depth is changed from... Unified conversion to The density of the sedimentary layer was reduced from Unified conversion to Substitute the values ​​into the calculation.

[0133] The lithospheric deflection and isostatic compensation sub-model is used to describe the effects of sedimentary load, crustal extension and thinning, and magma underplating on the basement interface, Moho discontinuity, and lithospheric interface burial depth. Under the local isostatic approximation, the vertical interface response caused by the load is expressed as:

[0134] ;

[0135] in, The vertical response of the interface, in units of ; The load per unit area is expressed in units of 1000 ppm. ; The density difference between the upper and lower media involved in the compensation is expressed in units of... .

[0136] When considering the elastic deflection of the lithosphere, the load response is expressed using the deflection equation of a thin elastic plate:

[0137] ;

[0138] in, Lithosphere flexural stiffness, unit: ; For mantle density, The density of seawater is given in units of 1 / 2. ; The vertical response of the interface, i.e., the flexural displacement, is expressed in units of 1000 ppm. .

[0139] The flexural stiffness Represented as:

[0140] ;

[0141] in, Young's modulus, unit: ; Effective elastic thickness, in units of ; Poisson's ratio is dimensionless. The above lithospheric flexure parameters are determined based on existing lithospheric flexure research results, passive continental margin basin flexure simulation results, or publicly available geodynamic parameter values ​​for the study area. Preferably, Young's modulus is used. Take 50-100 Effective elastic thickness Take 5000-20000 Poisson's ratio Take a value of 0.25-0.30, mantle density Take 3300-3350 Seawater density Take 1020-2030 .

[0142] It is important to note that magma underplating is represented in the model by a lower crustal high-velocity body material layer. This lower crustal high-velocity body has density, thermal conductivity, radioactive heat generation rate, and viscosity parameters that differ from those of ordinary lower crust. Its geometric distribution is defined by the top and bottom interfaces of the lower crustal high-velocity body, and its physical properties are determined based on magmatic activity data or existing rock property studies. The lower crustal high-velocity body influences the mass distribution and isostatic compensation results near the Mohorovičić discontinuity through its density difference, and affects the temperature field and heat flow distribution of the lower crust and lithosphere through its thermal property differences.

[0143] By coupling the crustal extension thinning sub-model, the lithospheric thermal evolution sub-model, the sedimentary load sub-model, and the deflection and isostatic compensation sub-model, a coupled numerical model of extension thinning-thermal evolution-interface burial depth response is formed. The Underworld2 numerical simulation platform is used, providing a unified two-dimensional mesh and particle material field. Specifically, the current crustal structure model provides the geometric framework and material property distribution; the crustal extension thinning sub-model provides the crustal thinning amount and extension factor; the lithospheric thermal evolution sub-model provides the temperature field and heat flow response; the sedimentary load sub-model provides the sedimentary loading amount; the deflection and isostatic compensation sub-model calculates the burial depth response of the basement interface, Moho discontinuity, and lithospheric interface based on load changes and density structure; and the magma underplating response expression is used to describe the influence of high-velocity anomalies at the bottom of the lower crust on density structure, thermal structure, and interface response.

[0144] The inputs to the coupled numerical model include current crustal geometric data, material properties, fault interpretation data, sedimentary layer thickness data, heat flow data, magmatic activity data, and Moho depth data. The outputs of the coupled numerical model include simulated crustal thickness, simulated basement depth, simulated Moho depth, and spatial distribution of high-velocity bodies in the lower crust. These outputs are used for simulations, objective function construction, and parameter updates in subsequent steps.

[0145] It should be noted that in this step, the initial crustal thickness... As input variables and variables to be updated in the coupled numerical model, they participate in the model expression, but their upper and lower bound constraints, different parameter ranges of different constructions, objective function construction and optimization iteration process are completed in subsequent steps.

[0146] Step (3): Based on the interpretation of seismic reflection profiles and drilling data, the depth of the same rift basement and the depth of the sedimentary bottom boundary after decompaction correction are used to establish the upper and lower limits of the initial crustal thickness.

[0147] The decompaction correction uses a holistic decompaction model to restore the sedimentary layers in the study area. This means that the sedimentary layers in the study area are considered as a whole, and the thickness of the sedimentary layers after decompaction is calculated based on the current total thickness, average porosity, and initial porosity, so as to eliminate the influence of sediment compaction on the interface burial depth and sediment load.

[0148] Specifically, let the current thickness of the sedimentary layer be... The thickness of the deposited layer after decompaction Represented as:

[0149] ;

[0150] in, and These represent the thickness of the sedimentary layer after decompaction and the current thickness of the sedimentary layer, respectively, in units of... ; and These are the average porosity and initial porosity of the sedimentary layer, respectively, and are dimensionless.

[0151] Furthermore, the average porosity is obtained statistically from drilling data, while the initial porosity is obtained by inversion using a compaction model with exponential decay. The compaction model is as follows:

[0152] ;

[0153] in, For burial depth Porosity at a given location, dimensionless; The unit is ; This is the compaction coefficient, in units of... e is the natural constant; generally, the compaction coefficient... Preferred concentration: 0.3-0.5 .

[0154] In this embodiment, a representative burial depth is selected. And by combining the current sedimentary layer thickness and average porosity, the initial porosity can be inverted, that is:

[0155]

[0156] To avoid excessive corrections in local thick depositional units from the overall decompaction model, physical constraints are imposed on the inversion parameters: the initial porosity is limited to a certain value. And the thickness of the deposited layer after decompaction meets the requirements. When the calculation result exceeds the above constraints, the result exceeding the limit will be truncated or the corresponding parameter combination will be removed.

[0157] Based on the seismic interpretation data and drilling data of the study area, the current sedimentary layer thickness and average porosity data of several typical geological units were statistically analyzed. The initial porosity and the sedimentary layer thickness after decompaction of different typical geological units in the study area were calculated using the above-mentioned overall decompaction model method. The results are shown in Table 3.

[0158] Table 3 Initial porosity and decompacted sedimentary layer thickness of typical geological units in the study area

[0159]

[0160] Based on the interpretation of seismic reflection profiles and the calculation of the predicted burial depth of key interfaces at each period using drilling data, the burial depth of the same rift basement is used as the lower limit constraint for the initial crustal thickness inversion, and the burial depth of the sedimentary bottom boundary after decompaction correction is used as the upper limit constraint for the initial crustal thickness inversion.

[0161] Based on the sedimentary thickness and density of each period, the burial depth of the syn-rift basement and the burial depth of the decompaction-corrected sedimentary bottom boundary are calculated. The lower limit of the burial depth of the syn-rift basement is used as the minimum true initial crustal thickness constraint, and the upper limit of the burial depth of the decompaction-corrected sedimentary bottom boundary is used as the maximum true initial crustal thickness constraint. This yields the upper and lower limits of the initial crustal thickness in the study area. In cases of significant magma underplating, the coupled numerical model is processed according to the magma intrusion layer as follows: if it is located at the bottom of the lower crust, near the Moho discontinuity, it is treated as a lower crustal high-velocity body; if it is located in the upper crust or sedimentary layer, it is treated as either the upper crust or sedimentary layer, respectively. This avoids misinterpreting the abnormal thickness of magma at different layers as the initial crustal thickness.

[0162] Step (4): Divide the study area into tectonic zones and set the initial crustal structure and geophysical parameter ranges for different tectonic zones;

[0163] Based on tectonic features, the study area is divided into tectonic zones, sequentially from continent to ocean: near-end zone, strongly thinning zone, far-end zone, ocean-continent transition zone (COT), and ocean basin. For each tectonic zone, and considering its crustal structure and tectonic characteristics, ranges are set for initial crustal thickness, current crustal thickness, thermolithospheric thickness, current upper crustal thickness, current lower crustal thickness, and the densities of sedimentary layers, upper crust, and lower crust. The near-end zone is represented by the Enping Depression, the strongly thinning zone by the Baiyun Depression (with the Baiyun East Depression selected as a typical modeling area), and the far-end zone by the Liwan Depression. The study area also includes the outer highlands extending towards the ocean. Based on the aforementioned upper and lower limits of initial crustal thickness, corresponding initial crustal thickness search ranges are set for each zone; for example, the near-end zone is set to 20-32. The strong reduction and thinning zone is set to 10-20. The remote band is set to 8-10. .

[0164] The relevant geophysical parameter ranges for different tectonic zone zoning units should be initialized according to Table 4. See Table 4 for details:

[0165] Table 4. Range of relevant geophysical parameters for tectonic zone zoning units.

[0166]

[0167] It should be noted that the parameter ranges listed in Table 4 above are used for subsequent forward simulation and iterative inversion of the coupled numerical model.

[0168] Step (5): Set the initial crustal thickness as the parameter to be inverted, perform forward modeling based on the coupled numerical model, and obtain the current simulation results;

[0169] The model, when performing forward modeling, considers the effects of fault displacement, lithospheric deflection, sedimentary load, and magma underplating on crustal structure. The forward modeling results are as follows: Figure 7 As shown.

[0170] Step (6): Construct the objective function based on the difference between the current simulation results and the measured constraint data;

[0171] The objective function is a weighted sum of at least two of the following: simulated crustal thickness error, simulated Moho discontinuity depth error, simulated heat flow error, and simulated basement burial depth error.

[0172] The initial crustal thickness was set as the parameter to be inverted. Forward simulation was performed based on each parameter combination of the constructed coupled numerical model of extensional thinning-thermal evolution-interface burial depth response to obtain the current simulated Moho depth, simulated basement depth, and simulated distribution of high-velocity bodies in the lower crust. It is important to note that while performing forward modeling, pre-existing faults were considered to characterize tectonic segmentation and differential extensional boundaries, lithospheric deflection reflects the elastic bending response caused by changes in tectonic load, sedimentary load is used to simulate the impact of sedimentary layer thickening on crustal subsidence, and magma underplating is used to characterize the impact of high-velocity bodies and thermal structural anomalies in the lower crust, thereby improving the consistency between the forward modeling results and current observational constraints.

[0173] An objective function is constructed based on the differences between the current simulation results and the measured constraint data. This objective function is a weighted sum of the errors in simulated crustal thickness and simulated basement depth. A hybrid optimization algorithm combining particle swarm optimization (PSO) and local gradient correction is used to iteratively solve for the inversion parameters. The PSO is used for global search, while the local gradient correction is used to fine-tune the parameters within the optimal region to improve the inversion convergence speed and solution accuracy. Specific formulas and explanations are as follows:

[0174] Let the parameter vector to be inverted be The objective function is The objective function is dimensionless and is a weighted sum of the simulated crustal thickness error and the simulated basement depth error.

[0175] ;

[0176] in, and These are the simulated crustal thickness and the measured crustal thickness, respectively, with units of... ; and These are the simulated and measured foundation depths, respectively, with units of [units not specified]. ; and , which is a weighting coefficient, dimensionless.

[0177] Let the first The particle in the first The velocity vector and position vector of the line are respectively and Its optimal position is The global optimal position is Since it is dimensionless, the update formula for the particle swarm optimization algorithm is:

[0178] ;

[0179] ;

[0180] in, For inertial weights, For individual learning factors, As a group learning factor, and It is a random number between 0 and 1, dimensionless.

[0181] When the particle swarm optimization algorithm iterates into the region where the objective function is relatively optimal, the parameters are further updated using the local gradient correction formula:

[0182] ;

[0183] in, and These represent the local optimal parameters at the th... generation and first The value of the substitute, Step size, The gradient of the objective function with respect to the parameters to be inverted is dimensionless.

[0184] During the iteration process, a global search is first performed using the particle swarm optimization algorithm. When the objective function decreases by less than a preset decrease threshold for several consecutive iterations, local gradient correction is initiated to refine the current optimal solution. This process is repeated until the change in the objective function is less than a preset convergence threshold or the preset maximum number of iterations is reached. The convergence condition is defined as the relative rate of change of the objective function in two adjacent iterations being less than a preset threshold, specifically expressed as:

[0185] ;

[0186] in, and They represent the first Second and third The objective function value of the next iteration. The threshold value is dimensionless.

[0187] When the above conditions are continuously less than the preset decreasing threshold When the above conditions are less than the preset convergence threshold, local gradient correction is initiated; when the above conditions are less than the preset convergence threshold, local gradient correction is initiated. Or reach the preset maximum number of iterations Stop iterating when the time is right. , All are dimensionless parameters and satisfy the following conditions: > > >0.

[0188] Specifically, a hybrid optimization algorithm combining particle swarm optimization and local gradient correction is used to invert the initial crustal thickness parameters. The relative rate of change of the objective function between two adjacent iterations is used as the iteration state criterion. Preferably, a preset threshold is used. Set the value to 0.01, and set the preset threshold for the decrease. Set the threshold to 0.005. Taking 0.0001 satisfies > > >0.

[0189] When the relative rate of change of the objective function is less than the preset decreasing threshold for 5 consecutive iterations When the particle swarm optimization algorithm is considered to have entered a locally optimal region where the objective function changes slowly, local gradient correction is initiated to further refine the current globally optimal parameters; when the relative rate of change of the objective function between two consecutive iterations is less than a preset convergence threshold... If the inversion result is considered to have reached the preset convergence accuracy, the iteration stops. If the number of iterations reaches the preset maximum number of iterations before the above convergence condition is met... Preferably, If the iteration is repeated 500 times, the iteration terminates, and the parameter combination corresponding to the minimum current objective function value is taken as the inversion result.

[0190] It is important to note that a manual interactive correction mechanism is set up during the parameter update process to receive correction instructions for abnormal parameters of local models, such as fracture location, layer boundary, density value or initial thickness range. After receiving the correction instruction, the system only performs local recalculation on the correction area and its adjacent areas, thereby improving the rationality of the recovery results under complex structural conditions, and continues to perform iterative inversion after correction.

[0191] Step (7): Update the parameters to be inverted using the optimization iterative algorithm, and repeat steps (5) to (6) until the preset convergence condition is met;

[0192] The optimization iterative algorithm includes particle swarm optimization, gradient descent, quasi-Newton algorithm, genetic algorithm, simulated annealing algorithm or a combination thereof; the parameters to be inverted include at least one of the following: initial crustal thickness zoning parameters, fault extension parameters, and magma underplating thickness parameters.

[0193] A manual interaction correction mechanism is set up to receive user correction instructions for local model parameters and continue iterative inversion after correction.

[0194] Step (8): Output the initial crustal thickness recovery results of the northern continental margin of the South China Sea.

[0195] The output includes an initial crustal thickness distribution map and its confidence interval for the northern continental margin of the South China Sea, as well as a tectonic zone extension coefficient distribution map. The confidence interval for the initial crustal thickness characterizes the uncertainty of the inversion results. It is calculated by retaining the feasible solution set obtained from multiple inversions, provided that preset convergence conditions are met, statistically analyzing the initial crustal thickness results for each region, and taking its quantile interval as the confidence interval. The tectonic zone extension coefficient distribution map characterizes the crustal extension intensity of tectonic zones in different regions. Specifically, the confidence interval is expressed as follows:

[0196] ;

[0197] in, For position The confidence interval for the initial crustal thickness. and These represent the lower and upper quantiles of the initial crustal thickness inversion result at this location, respectively. For the preset significance level, preferably, The value is usually taken as 0.05.

[0198] Specifically, its elongation coefficient is expressed as:

[0199] ;

[0200] in, is the extension coefficient of the tectonic zone, dimensionless; This represents the initial crustal thickness. The current thickness of the Earth's crust, in units of 1 / 2. Based on the initial crustal thickness obtained from the inversion of each tectonic zone and the current crustal thickness, the corresponding extension coefficient is calculated, and a distribution map of the tectonic zone extension coefficient is generated.

[0201] Taking the Baiyun-Liwan Depression as an example, The value is 0.05. Output the initial crustal thickness reconstruction results of the Baiyun-Liwan Depression on the northern continental margin of the South China Sea and its ocean-facing outer highlands, including the initial crustal structure profile of the Baiyun-Liwan Depression reconstructed based on the Underworld2 coupled numerical model, as shown below. Figure 8 As shown, and the initial crustal thickness and confidence interval (95% CI) and extension coefficient β distribution maps of the Baiyun-Liwan Depression, which were inverted and restored based on the Underworld2 coupled numerical model, are as follows. Figure 9 As shown in the figure, the inversion results indicate that the initial crustal thickness in the study area decreases from land to ocean. This trend is consistent with the thinning pattern of the northern passive continental margin of the ancient South China Sea before the start of the new South China Sea rift, and also confirms the general view that the lithosphere had already undergone a certain degree of thinning before the rift. Meanwhile, the distribution results of the tectonic zone extension coefficient β further show that the restored values ​​of different tectonic zones differ significantly, with greater extension corresponding to the present-day ocean-continent transition zone (COT) and the ocean-bound basin area, while the extension of the relatively land-side area is smaller. Furthermore, the simulations of the current crustal thickness, Moho depth, and distribution of the lower crustal high-velocity body by the Underworld2 coupled numerical model are in good agreement with the current observed values ​​of the northern South China Sea continental margin, indicating that the method can effectively restore the initial crustal thickness of the northern South China Sea continental margin, and the restoration results are reasonable. Compared with existing technologies, this invention improves the accuracy and stability of the initial crustal thickness recovery results by unifying multi-source geophysical and geological data under the same inversion framework and combining a hybrid optimization strategy of particle swarm optimization (PSO) and local gradient correction. Furthermore, it quantitatively characterizes the inversion uncertainty through a 95% confidence interval (95% CI), thereby enhancing the scientific validity and reliability of the results.

Claims

1. A method for restoring the initial crustal thickness of a passive continental margin, characterized in that, Includes the following steps: Step (1): Obtain multi-source geophysical and geological data, including seismic reflection profile interpretation data, present crustal thickness data, Moho discontinuity depth data, drilling data, magma data, heat flow data, sedimentary layer thickness data, basement depth data, and tectonic zone delineation data. Step (2): Based on the multi-source geophysical and geological data, establish a current crustal structure model and a coupled numerical model of extensional thinning-thermal evolution-interface burial depth response; Step (3): Based on the interpretation of seismic reflection profiles and drilling data, the depth of the same rift basement and the depth of the sedimentary bottom boundary after decompaction correction are used to establish the upper and lower limits of the initial crustal thickness. Step (4): Divide the study area into tectonic zones and set the initial crustal structure and geophysical parameter ranges for different tectonic zones; Divide the study area into tectonic zones according to tectonic characteristics, from land to ocean: near-end zone, strong thinning zone, far-end zone, ocean-continent transition zone (COT), ocean basin; For the crustal structure and tectonic characteristics of different tectonic zones, set the initial crustal thickness, current crustal thickness, thermolithospheric thickness, current upper crustal thickness, current lower crustal thickness, and the range of values ​​for sedimentary layer, upper crust, and lower crustal density respectively; Among them, the near-end zone is represented by Enping Depression, the strong thinning zone is represented by Baiyun Depression and Baiyun East Depression is selected as a typical modeling area, the far-end zone is represented by Liwan Depression, and the study area also includes the outer edge highland area in the direction of ocean; Step (5): Set the initial crustal thickness as the parameter to be inverted, and perform forward simulation based on the coupled numerical model to obtain the current simulation results; the model considers the influence of fault displacement, lithosphere deflection, sedimentary load and magma underplating on the crustal structure while performing forward simulation. Step (6): Construct an objective function based on the difference between the current simulation results and the measured constraint data; the objective function is the weighted sum of at least two of the following: simulated crustal thickness error, simulated Moho depth error, simulated heat flow error, and simulated basement depth error. Step (7): Update the parameters to be inverted using an optimization iterative algorithm, and repeat steps (5) to (6) until the preset convergence condition is met; the optimization iterative algorithm includes particle swarm optimization, gradient descent, quasi-Newton algorithm, genetic algorithm or simulated annealing algorithm, or a combination of any two of the iterative algorithms; the parameters to be inverted include at least one of the initial crustal thickness zoning parameters, fault extension parameters or magma underplating thickness parameters; Step (8): Output the initial crustal thickness restoration results of the northern continental margin of the South China Sea: The specific steps in step (2) above for constructing the coupled numerical model of stretching thinning, thermal evolution, and interface burial depth response are as follows: Based on the existing crustal structure model, a coupled numerical model of extensional thinning, thermal evolution, and interface burial depth response is further established. This coupled numerical model includes a crustal extensional thinning sub-model, a lithospheric thermal evolution sub-model, a sedimentary load sub-model, and a deflection and isostatic compensation sub-model. By coupling these sub-models, a coupled numerical model of extensional thinning, thermal evolution, sedimentary load, and deflection and isostatic compensation is formed. Specifically, the existing crustal structure model provides the geometric framework and material property distribution; the crustal extensional thinning sub-model provides the crustal thinning amount and extension factor; the lithospheric thermal evolution sub-model provides the temperature field and heat flow response; the sedimentary load sub-model provides the sedimentary loading amount; the deflection and isostatic compensation sub-model calculates the burial depth response of the basement interface, Moho discontinuity, and lithospheric interface based on load changes and density structure; and the magma underplating response expression is used to describe the influence of high-velocity anomalies at the bottom of the lower crust on density structure, thermal structure, and interface response. The inputs to the coupled numerical model include current crustal geometric data, material property parameters, fault interpretation data, sedimentary layer thickness data, heat flow data, magmatic activity data, and Moho depth data. The output of the coupled numerical model includes simulated crustal thickness, simulated basement depth, simulated Moho depth, and spatial distribution of high-velocity bodies in the lower crust.

2. The method for restoring the initial crustal thickness of a passive continental margin according to claim 1, characterized in that: The specific steps for constructing the current crustal structure model in step (2) are as follows: First, based on step (1), according to the tectonic trend of the study area, the distribution direction of the seismic reflection profile, and the research objectives, a two-dimensional research profile passing through the Baiyun East Depression, the Liwan Depression, and the ocean-facing outer high region of the study area is selected as the modeling profile; a two-dimensional rectangular coordinate system is established along the modeling profile, wherein the horizontal direction is... Direction, vertical depth direction is Direction, among which Corresponding to mean sea level, The area is the air layer. The area consists of the seawater layer and the underlying geological bodies; Based on seismic reflection profile interpretation data, well layer data, current crustal thickness data, Moho depth data, and basement depth data, the main geological interface control points on the modeling profile are extracted. The main geological interfaces include at least the sea level, seafloor interface, sedimentary basement interface, Conrad discontinuity, Moho discontinuity, and lithospheric boundary. The current crustal structure model is divided into eight layers from top to bottom: air layer, seawater layer, sedimentary layer, upper crust, lower crust, lower crustal high-velocity body, lithospheric mantle, and asthenosphere.

3. The method for restoring the initial crustal thickness of a passive continental margin according to claim 1, characterized in that: Specifically, in step (3): let the current thickness of the sedimentary layer be... The thickness of the deposited layer after decompaction Represented as: ； in, and These represent the thickness of the sedimentary layer after decompaction and the current thickness of the sedimentary layer, respectively, in units of... ; and These are the average porosity and initial porosity of the sedimentary layer, respectively, and are dimensionless. Furthermore, the average porosity is obtained statistically from drilling data, while the initial porosity is obtained by inversion using an exponentially decaying compaction model. The compaction model is as follows: ； in, For burial depth Porosity at a given location, dimensionless; The unit is ; This is the compaction coefficient, in units of... e is the natural constant; compaction coefficient It is 0.3-0.5 ; Based on the interpretation of seismic reflection profiles and the calculation of the predicted burial depth of key interfaces at each period using drilling data, the burial depth of the same rift basement is used as the lower limit constraint for the initial crustal thickness inversion, and the burial depth of the sedimentary bottom boundary after decompaction correction is used as the upper limit constraint for the initial crustal thickness inversion.

4. The method for restoring the initial crustal thickness of a passive continental margin according to claim 1, characterized in that: Specifically in step (6): Let the parameter vector to be inverted be The objective function is The objective function is dimensionless and is a weighted sum of the simulated crustal thickness error and the simulated basement depth error. ； in, and These are the simulated crustal thickness and the measured crustal thickness, respectively, with units of... ; and These are the simulated and measured foundation depths, respectively, with units of [units not specified]. ; and These are weighting coefficients, dimensionless. Let the first The particle in the first The velocity vector and position vector of the line are respectively and Its optimal position is The global optimal position is Since it is dimensionless, the update formula for the particle swarm optimization algorithm is: ； ； in, For inertial weights, For individual learning factors, As a group learning factor, and A random number between 0 and 1, dimensionless; When the particle swarm optimization algorithm iterates into the region where the objective function is relatively optimal, the parameters are further updated using the local gradient correction formula: ； in, and These represent the local optimal parameters at the th... generation and first The value of the substitute, Step size, The gradient of the objective function with respect to the parameters to be inverted is dimensionless. During the iteration process, a global search is first performed using the particle swarm optimization algorithm. When the objective function decreases by less than a preset decrease threshold for several consecutive iterations, local gradient correction is initiated to refine the current optimal solution. This process is repeated until the change in the objective function is less than a preset convergence threshold or the preset maximum number of iterations is reached. The convergence condition is defined as the relative rate of change of the objective function in two adjacent iterations being less than a preset threshold, specifically expressed as follows: ； in, and They represent the first Second and third The objective function value of the next iteration. The threshold value is dimensionless. When the above conditions are continuously less than the preset decreasing threshold When the above conditions are less than the preset convergence threshold, local gradient correction is initiated; when the above conditions are less than the preset convergence threshold, local gradient correction is initiated. Or reach the preset maximum number of iterations Stop iterating when the time comes; , All are dimensionless parameters and satisfy the following conditions: > > >

0.

5. The method for restoring the initial crustal thickness of a passive continental margin according to claim 1, characterized in that: The output results in step (8) include the initial crustal thickness distribution and its confidence interval in the northern continental margin of the South China Sea, as well as the distribution map of tectonic zone extension coefficients; wherein, the confidence interval of the initial crustal thickness is used to characterize the uncertainty of the inversion results, and its calculation method is as follows: under the premise of satisfying the preset convergence conditions, retain the feasible solution set obtained by multiple inversions, statistically analyze the initial crustal thickness results of each region, and take its quantile interval as the confidence interval; while the distribution map of tectonic zone extension coefficients is used to characterize the crustal extension intensity of tectonic zones in different regions; specifically, the confidence interval is expressed as: ； in, For position The confidence interval for the initial crustal thickness. and These represent the lower and upper quantiles of the initial crustal thickness inversion result at this location, respectively. At the preset significance level, Take 0.05; Specifically, its elongation coefficient is expressed as: ； in, is the extension coefficient of the tectonic zone, dimensionless; The initial crustal thickness, The current thickness of the Earth's crust, in units of 1 / 2. .

Images