Method for inversion of direct current data using multi-point geostatistical structural constraints

Abstract

The application discloses a method for DC data inversion constrained by multi-point geostatistics structure, and belongs to the field of geophysical exploration, which comprises the following steps: based on geological drilling data, a model containing geological structure information is established by using multi-point geostatistics algorithm, and an underground stratum structure operator is constructed by using the model; the structure operator is added to a regularization term, a data residual term is calculated by giving an initial resistivity model, and an inversion objective function with structure constraint is formed; in the inversion process, a limited memory quasi-Newton optimization algorithm is used to iteratively optimize the objective function until the difference between the forward simulation data and the measured data reaches a set threshold, that is, the constrained inversion is completed; and finally, the resistivity inversion model is output. The application reduces the problems of stratum boundary blurring and resolution decline caused by the regularization smoothing constraint in the traditional method, is suitable for underground stratum structure detection and fine imaging under complex geological environment, and has wide application value.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of geophysical exploration, in particular to a method for constraining direct current method data inversion by using multi-point geostatistical structure. BACKGROUND

[0002] With the rapid advancement of urban construction, the accuracy requirements for underground space exploration are continuously improved, especially in the fields of urban infrastructure construction, active fault monitoring, underground resource exploration, etc., the demand for high-resolution, non-destructive exploration technology is growing. Direct current method is a geophysical exploration method for underground structure, which collects direct current method data and constructs resistivity model by inverting direct current method data. Direct current method technology has significant advantages in urban environment with flexible observation mode, high resolution and strong anti-interference ability. In practical application, the inversion problem is usually ill-posed, that is, the observation data is limited and contains noise, which leads to instability of the inversion solution. The existing technology mainly relies on mathematical regularization to optimize the solution, but these technologies only constrain the smoothness of the solution from the numerical point of view, and fail to fully utilize the prior structural information of the underground medium contained in the drilling data, resulting in low resolution of the inversion result. SUMMARY

[0003] The present application aims to overcome the shortcomings of the prior art and provide a method for constraining direct current method data inversion by using multi-point geostatistical structure, which introduces prior structural constraints in the inversion process, so that the inversion solution meets the geological structure characteristics while being optimized mathematically, improving the resolution and stability of imaging.

[0004] The present application aims to overcome the shortcomings of the prior art and provide a method for constraining direct current method data inversion by using multi-point geostatistical structure, which introduces prior structural constraints in the inversion process, so that the inversion solution meets the geological structure characteristics while being optimized mathematically, improving the resolution and stability of imaging.

[0005] Underground formation structure operator construction: based on geological drilling data, multi-point geostatistical algorithm is used to randomly model the underground space structure, and a multi-point geostatistical model containing geological structure information is generated, and an underground formation structure operator is constructed based on the multi-point geostatistical model;

[0006] Establishment of structure-constrained inversion objective function: the structure operator is added to the regularization term, and the data residual term is calculated based on the initial resistivity model to form the inversion objective function with structure constraint;

[0007] Structure-constrained inversion: the finite memory quasi-Newton optimization algorithm is used to iteratively optimize the inversion objective function with structure constraint until the difference between the forward simulation data and the measured data reaches the preset constraint condition, and the inversion model is obtained;

[0008] Inversion model output: the inversion model is output for geological interpretation.

[0009] Further, the multi-point geostatistical algorithm is a direct sampling algorithm, which adopts a flexible and variable search neighborhood to determine data events, is suitable for non-uniform media, and has high calculation efficiency.

[0010] Further, the constructing an underground formation structure operator based on the multi-point geostatistical model specifically includes: performing grid division on the multi-point geostatistical model, calculating information entropy of a lithology or physical property parameter corresponding to each grid, and constructing a normalized structure operator based on the information entropy.

[0011] Further, the constructing an underground formation structure operator based on the multi-point geostatistical model further includes: constructing a structure operator matrix based on the structure operator, dynamically weighting a smoothing degree of the model according to a distribution of information entropy of the lithology or physical property parameter in the multi-point geostatistical model, and reducing a weight of a smoothing constraint in a region with high information entropy value, the region with high information entropy value including a geological interface and a structure mutation region.

[0012] Further, the regularization term is a difference between a resistivity inversion model and a reference model, and is constructed by using an L2 norm, wherein the reference model is a uniform model with an average resistivity of a detection region.

[0013] Further, the data residual term is a difference between measured data and forward simulation data, and is constructed by using an L2 norm, wherein the measured data is an actual measurement value, and the forward simulation data is a theoretical prediction value corresponding to a current model.

[0014] Further, the inversion objective function with a structure constraint includes a data residual term and a structure constraint term, and a regularization term is added to a weight of the inversion objective function to prevent overfitting.

[0015] Further, the structure constraint inversion process adopts a multi-iteration strategy to update the model, and a termination condition is set as: a norm formed by a difference between the measured data and the forward simulation data multiplied by a data weighting matrix, when a value of the termination condition is less than a set threshold or a change rate meets a convergence condition, the iteration is stopped; wherein the data weighting matrix is used to represent importance of the observation value.

[0016] Further, the structure constraint inversion process includes: updating a model parameter based on information of a gradient of the objective function, and setting a parameter conversion relationship between the model parameter and an actual physical property; the conversion relationship is a nonlinear function, which is used to map the parameter obtained by the inversion to a corresponding physical property value, and further constrain a range of the physical property parameter.

[0017] The application further provides an electronic device comprising a memory and a processor, wherein the memory is coupled to the processor; the memory is used for storing program data, and the processor is used for executing the program data to realize the method.

[0018] The application further provides a computer readable storage medium, which stores a computer program, and the program is executed by a processor to realize the method.

[0019] The application has the beneficial effect that the application combines the prior information of the formation structure provided by drilling and the like, establishes a model containing the geological structure information based on the multiple-point geostatistics, and uses the model to construct the underground formation structure operator; the formation structure operator is added to the regularization term to effectively constrain the formation structure in the direct current method data inversion process, and the resolution of the resistivity inversion model is improved. The above technical means effectively improves the resolution of the direct current method data inversion result, and especially under the condition of low signal-to-noise ratio or uneven data coverage, still has good convergence and imaging quality. Reduces the problems of blurred formation boundary and reduced resolution caused by the regularization smoothing constraint of the traditional method. It is suitable for urban underground structure detection and fine imaging in complex geological environment, and has wide application value. BRIEF DESCRIPTION OF DRAWINGS

[0020] In order to more clearly illustrate the technical solutions in the embodiments of the application, the following will briefly introduce the drawings needed to be used in the embodiment description. Obviously, the drawings in the following description are only some embodiments of the application, and other drawings can be obtained by those skilled in the art without creative labor.

[0021] Figure 1 is a flow chart of the direct current method data inversion constrained by the multiple-point geostatistics structure;

[0022] Figure 2 is a schematic diagram of actual case inversion result, wherein the upper drawing is a schematic diagram of structure-constrained inversion result; the lower drawing is a schematic diagram of structure-constrained inversion result, and the red soil layer is shown as a light blue area in the drawing. DETAILED DESCRIPTION

[0023] The application will be described in detail below with reference to the drawings. In the case of no conflict, the features in the following embodiments and implementation manners can be combined with each other.

[0024] The embodiment of the application provides a method for direct current method data inversion constrained by the multiple-point geostatistics structure, as shown in Figure 1 The method comprises the following steps:

[0025] The underground formation structure operator is constructed: based on the geological drilling data, a multi-point geostatistical algorithm is used to randomly model the underground space structure, a multi-point geostatistical model containing geological structure information is generated, and the underground formation structure operator is constructed based on the multi-point geostatistical model.

[0026] The establishment of the structure-constrained inversion objective function: the structure operator is added to the regularization term, the data residual term is calculated by taking the uniform resistivity model with a resistivity value of 30 ohm meters as the initial model, and the inversion objective function with structure constraint is formed;

[0027] Structure-constrained inversion: the finite memory quasi-Newton optimization algorithm is used to iteratively optimize the inversion objective function with structure constraint until the difference between the forward simulation data and the measured data reaches the preset constraint condition, and the inversion model is obtained;

[0028] Inversion model output: the above inversion model is output and used for geological interpretation.

[0029] The multi-point geostatistical algorithm is a direct sampling algorithm, which uses a flexible and variable search neighborhood to determine data events, is suitable for non-uniform media, and has high computational efficiency.

[0030] In one embodiment, the formation structure operator is structure information quantitatively extracted from the multi-point geostatistical model using information entropy, specifically: assuming that the multi-point geostatistical model contains N kinds of lithology or physical properties, the model is grid divided, and the probability of the lithology or physical property corresponding to the i-th grid of the multi-point geostatistical model is represented by The information entropy of the lithology or physical property corresponding to the grid is :

[0031] (1);

[0032] The normalized structure operator is constructed using information entropy :

[0033] (2);

[0034] The structure operator matrix is established based on the structure operator , According to the distribution of the structure parameters, the smoothness of the resistivity inversion model is dynamically weighted, the weight of the smoothness constraint is reduced in the area with high information entropy value, and the area with high information entropy value includes geological interfaces and structure mutation areas, etc., so as to obtain clear structure information.

[0035] In one embodiment, the regularization term with the structure operator is in the form of:​

[0036] (3);

[0037] In formula (3), m is the resistivity inversion model, is the reference model, is the smoothing matrix, is the structure operator matrix.

[0038] The data residual term of the constraint inversion objective function is defined as , which describes the difference between the actual data and the model theoretical data:

[0039] (4);

[0040] In formula (4), m is the resistivity inversion model, is the theoretical data generated by forward simulation, is the actual data measured. For the direct current method, m and ( ) correspond to the model resistivity and apparent resistivity data, respectively. is the inverse of the actual data covariance matrix, which is used to reduce the influence of large error data.

[0041] Comprehensive formula (3) and formula (4), define the multi-point geostatistical structure constraint inversion objective function as,

[0042] (5);

[0043] In formula (5), , respectively represent the direct current method data residual term and the multi-point geostatistical structure constraint term of the resistivity inversion model. When performing constraint inversion, the structure constraint term in formula (5) can be weighted to adjust its weight.

[0044] The finite memory quasi-Newton method is used to solve the constraint inversion objective function (formula (5)). First, calculate the gradient of the constraint inversion objective function with respect to the resistivity inversion model, then apply the finite memory quasi-Newton method to convert the gradient into a search direction and an optimization step , and finally update the current resistivity inversion model for iteration, and the iteration calculation formula is,

[0045] (6).

[0046] In one embodiment, the inversion process includes using a multiple iteration update strategy and setting the following optimization convergence conditions:

[0047] .

[0048] in, Represents a data weighting matrix. Represents observation data, This represents forward simulation data. This indicates the number of data points.

[0049] In one embodiment, the inversion process includes parameter updates based on the gradient of the objective function, introducing a parameter transformation between the inversion parameters and the physical property (resistivity), with the transformation function as follows:

[0050] ;

[0051] , These are the minimum and maximum resistivity, respectively, used to ensure that any element of the inversion model vector... When the resistivity of the forward model changes, Constrained and Between. Because the resistivity values ​​vary over a large range, logarithmic coordinates were used to convert the resistivity parameters.

[0052] In one embodiment, the following specific experiments were conducted at the Liangzhu Ancient City Site in Hangzhou, Zhejiang Province, based on this method, including the following steps:

[0053] 1. Construction of Subsurface Stratigraphic Structure Operator: Based on geological drilling data, a model containing geological structural information is established using a multi-point geostatistical algorithm. This model is then used to construct the subsurface stratigraphic structure operator.

[0054] In this embodiment, the work area is located at the Liangzhu Ancient City Site in Hangzhou, Zhejiang Province. Drilling revealed a layer of burnt red soil approximately 0.2m thick and buried at a depth of about 2m, dating back to the Liangzhu period and presumably remnants of ancient architecture. The DC resistivity tomography (DCS) data acquisition electrode spacing was 0.5m, the survey line length was 24m, and data acquisition was performed using a Wenner device. To add structural constraints during the DCS data inversion process, two wells with a depth of 4m were drilled. A multi-point geostatistical method was used to calculate a model containing geological structures, and a stratigraphic structure operator was constructed using information entropy theory.

[0055] 2. Establishment of the structural constraint inversion objective function: Add the structural operator to the regularization term, calculate the residual term of the data given the initial model, and form an inversion objective function with structural constraints.

[0056] In this embodiment, the initial model adopted is a uniform model with a resistivity of 30 ohm-meters, and the initial model has a horizontal length of 24 meters and a vertical length of 4 meters.

[0057] 3. Structural Constraint Inversion: During the inversion process, the objective function is optimized by introducing structural constraint matrix operators into the inversion framework. These operators work together with traditional data terms and regularization terms to iteratively optimize the inversion process. The inversion process uses a finite-memory quasi-Newton method to iteratively optimize the objective function. A threshold value of 1 is set for the difference between the forward data and the measured data. When the difference between the forward data and the measured data reaches the set threshold value, the constraint inversion is completed.

[0058] In this embodiment, based on the acquisition method and drilling data, the inversion model size is selected as 24m × 4m, and the basic inversion grid is 0.1m × 0.1m. The initial inversion model adopts a uniform model with a uniformly distributed resistivity value of 30 ohm-meters.

[0059] The inversion process first uses the L-Curve method to determine the weights of the regularization term, with a final weight value of 100. Then, this weight is used as the weight of the regularization term in the structural constraint inversion objective function, and optimization inversion is performed simultaneously on both the structural constraint term and the data term. The inversion RMS error is then calculated using the system.

[0060] like Figure 2 As shown, the upper figure is the resistivity profile obtained by traditional inversion, and the lower figure is the resistivity profile corresponding to multi-point geostatistical structure-constrained inversion. From Figure 2 Comparing the results of the two inversion techniques, the thin layer of red-burnt soil was not displayed on the traditional inversion profile, while the multi-point geostatistical structure-constrained joint inversion technique was able to distinguish the thin layer of about 0.2m.

[0061] In summary, the DC electrical resistivity tomography (DCIP) data inversion method constrained by multi-point geostatistical structure can reduce the impact of traditional inversion regularization and smoothing constraints, and significantly improve the resolution of the inversion results while maintaining inversion accuracy. Compared with traditional inversion methods, the multi-point geostatistical structure constrained inversion technique is superior in terms of morphological characterization of subsurface thin-layer anomalies, accuracy of inverted physical property parameters, and accuracy of inverted thickness.

[0062] This invention also provides an electronic device, including a memory and a processor, wherein the memory is coupled to the processor; wherein the memory is used to store program data, and the processor is used to execute the program data to implement the method for inverting DC electrical resistivity tomography data constrained by multi-point geostatistical structures as described in any embodiment.

[0063] This invention also provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the method for inverting DC electrical resistivity tomography data constrained by multi-point geostatistical structures as described in any embodiment.

[0064] The computer readable storage medium can be an internal storage unit of any of the aforementioned data processing capable devices, such as a hard disk or a memory. The computer readable storage medium can also be any of the aforementioned data processing capable devices, such as a plug-in hard disk, a Smart Media Card (SMC), an SD card, a Flash Card, etc. Further, the computer readable storage medium can also include both an internal storage unit of any of the aforementioned data processing capable devices and an external storage device. The computer readable storage medium is used to store the computer program and other programs and data required by the aforementioned data processing capable devices, and can also be used to temporarily store data that has been output or will be output.

[0065] Those skilled in the art can understand that the above description is only preferred examples of the application and is not used to limit the application, although the application is described in detail with reference to the foregoing examples, those skilled in the art can still modify the technical solutions recorded in the foregoing examples or make equivalent replacements for part of the technical features. Any modification, equivalent replacement, etc. made within the spirit and principles of the application shall be included in the protection scope of the application.

Claims

1. A method for inverting DC resistivity tomography data using multi-point geostatistical structure constraints, characterized in that, Includes the following steps: Construction of Subsurface Stratigraphic Structure Operator: Based on geological drilling data, a multi-point geostatistical algorithm is used to stochastically model the subsurface spatial structure, generating a multi-point geostatistical model containing geological structural information. Subsurface stratigraphic structure operators are then constructed based on this multi-point geostatistical model. The stratigraphic structure operator is the structural information quantitatively extracted from the multi-point geostatistical model using information entropy. Specifically, assuming the multi-point geostatistical model contains N types of lithology or physical properties, the model is meshed, and... Let represent the probability of the lithology or physical property corresponding to the i-th grid in a multi-point geostatistical model, then the information entropy of the lithology or physical property corresponding to that grid. for: ； Using information entropy, a normalized structure operator is constructed. : ； Based on structure operators Establish the structure operator matrix , The smoothness of the resistivity inversion model will be dynamically weighted according to the distribution of structural parameters. In regions with high information entropy, the weight of the smoothness constraint will be reduced to obtain clear structural information. Establishment of the structurally constrained inversion objective function: The structural operator is added to the regularization term, and the residual term of the initial resistivity model calculation data is given to form the inversion objective function with structural constraints; specifically: Regularization terms with structured operators The format is: ; in, As a reference model, For smoothing matrices, The structure operator matrix; Define the data residuals of the constraint inversion objective function This is used to describe the difference between actual data and theoretical data from the model: ； Where m is the resistivity inversion model, These are theoretical data generated by forward modeling. These are actual data obtained through measurement; for the DC resistivity method, m is the corresponding model resistivity. and This corresponds to the apparent resistivity data of the model; It is the inverse of the covariance matrix of the actual data, used to reduce the impact of large error data; Define the objective function for multi-point geostatistical structure constraint inversion. for: ; in, , These represent the residual terms of the DC resistivity method data and the multi-point geostatistical structure constraint terms of the resistivity inversion model, respectively. During constraint inversion, the structural constraint terms in the inversion objective function can be weighted to adjust their weights. Structural constraint inversion: The inversion objective function with structural constraints is iteratively optimized using a finite-memory quasi-Newton optimization algorithm until the difference between the forward simulation data and the measured data reaches the preset constraint condition, thus obtaining the inversion model; Inversion model output: The inversion model is output for geological interpretation.

2. The method according to claim 1, characterized in that, The specific steps for constructing the underground stratigraphic structure operator based on the multi-point geostatistical model are as follows: the multi-point geostatistical model is divided into grids, the information entropy of the lithology or physical property parameters corresponding to each grid is calculated, and a normalized structure operator is constructed based on the information entropy.

3. The method according to claim 2, characterized in that, The method of constructing the underground stratigraphic structure operator based on the multi-point geostatistical model further includes: constructing a structure operator matrix based on the structure operator; dynamically weighting the smoothness of the model according to the distribution of information entropy of lithological or physical property parameters in the multi-point geostatistical model; reducing the weight of smoothing constraints in regions with high information entropy values; and the regions with high information entropy values ​​include geological interfaces and structural abrupt change regions.

4. The method according to claim 1, characterized in that, The regularization term represents the difference between the resistivity inversion model and the reference model, and is constructed using the L2 norm. The reference model is a uniform model whose resistivity value is the average resistivity of the probe region.

5. The method according to claim 4, characterized in that, The data residual term represents the difference between the measured data and the forward simulation data, constructed using the L2 norm. The measured data refers to the actual measured values, while the forward simulation data refers to the theoretical predicted values ​​corresponding to the current model.

6. The method according to claim 1, characterized in that, The inversion objective function with structural constraints includes a data residual term and a structural constraint term. Overfitting is prevented by increasing the weight of the regularization term in the inversion objective function.

7. The method according to claim 1, characterized in that, The structural constraint inversion process employs a multi-iteration strategy for model updates. The termination condition is set as follows: the norm formed by the difference between the measured data and the forward simulation data is multiplied by the data weighting matrix. When its value is less than a set threshold or the rate of change meets the convergence condition, the iteration stops. The data weighting matrix is ​​used to characterize the importance of the observed values.

8. The method according to claim 1, characterized in that, The structural constraint inversion process includes: updating the model parameters based on the information of the objective function gradient, and setting a parameter transformation relationship between the model parameters and the actual physical properties; the transformation relationship is a nonlinear function used to map the inverted parameters to the corresponding physical property values, and further constrain the range of physical property parameters.

9. An electronic device comprising a memory and a processor, characterized in that, The memory is coupled to the processor; wherein the memory is used to store program data, and the processor is used to execute the program data to implement the method as described in any one of claims 1-8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method as described in any one of claims 1-8.

Images