A time-lapse inversion method based on time-lapse term constraints

Abstract

The present application provides a time-lapse inversion method based on time-lapse term constraints, including: obtaining apparent resistivity data; dividing the research area of the underground medium into multiple grid cells, the initial resistivity value of each grid cell being determined according to an initial resistivity model, constructing an objective function, the objective function including a time-lapse term and a first weight factor, the first weight factor being the product of a first weight coefficient and a second weight coefficient, the first weight coefficient being used to reflect the resistivity change amount of the resistivity model at adjacent times at each grid cell, and the second weight coefficient being used to reflect the change amount of the resistivity value of the resistivity model at the current time at each grid cell relative to the resistivity value at an adjacent grid cell; using the apparent resistivity data, continuously iterating the model parameters of the resistivity model to continuously reduce the objective function, and obtaining an optimal resistivity model. The present application selects an appropriate weight factor for the time-lapse term to reduce the error of the inversion result and improve the accuracy of the inversion result.

Description

Technical Field

[0001] This application relates to the technical field of geophysical exploration. Specifically, it relates to a time-lapse inversion method based on time-lapse term constraint. Background Art

[0002] Geophysics is one of the main disciplines of earth science and is a comprehensive discipline for studying the earth and searching for mineral resources inside the earth. According to the different physical fields studied, it can be divided into gravity exploration, magnetic exploration, electrical exploration, seismic exploration, etc. The resistivity method is one of the main branch methods of electrical exploration. Due to its advantages of simple operation, low cost, high working efficiency, etc., it is widely used in scenarios such as environmental protection, ecological restoration, seawater intrusion, and geological disasters.

[0003] Since the changes in the external environment during the monitoring process affect the monitoring data, resulting in the inversion results of the monitoring data not being able to accurately invert the changes in the underground medium, the time-lapse inversion method based on time-lapse term constraint is often used in the prior art for inversion. Among them, the time-lapse term synchronously inverts the data at different times, and the data at different times constrain each other, which can ensure that the inversion results at different times change continuously with time and has a certain effect of suppressing background noise interference. However, in previous studies, the weight factor of the time-lapse term was usually set as a constant. This method ignores the situation that some areas underground have changed greatly within a certain period of time rather than changing continuously with time, which easily leads to a certain error in the inversion results and makes the accuracy of the inversion results relatively low. Summary of the Invention

[0004] In view of this, the purpose of this application is to provide a time-lapse inversion method based on time-lapse term constraint, which reduces the error of the inversion results by selecting a suitable weight factor for the time-lapse term, thereby improving the accuracy of the inversion results.

[0005] In a first aspect, an embodiment of this application provides a time-lapse inversion method based on time-lapse term constraint. The time-lapse inversion method includes:

[0006] Obtain the apparent resistivity data of the underground medium;

[0007] Divide the research area corresponding to the underground medium into multiple grid cells, and the initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model. The spatial range of the research area is larger than the spatial range of the observation area of the underground medium;

[0008] Construct an objective function for time-lapse inversion, where the objective function includes a time-lapse term and a first weight factor corresponding to the time-lapse term. The first weight factor is used to adjust the proportion of the time-lapse term in the objective function. Among them, the first weight factor is the product of a first weight coefficient and a second weight coefficient. The first weight coefficient is used to reflect the change in resistivity at each grid cell of the resistivity model at adjacent times, and the second weight coefficient is used to reflect the change in the resistivity value of the resistivity model at each grid cell at the current time relative to the resistivity value at adjacent grid cells;

[0009] Utilize the apparent resistivity data to continuously iteratively update the model parameters of the resistivity model to continuously reduce the objective function, so as to obtain an optimal resistivity model.

[0010] In an optional embodiment, the first weight coefficient is determined through the following steps:

[0011] Construct a resistivity change matrix, where the matrix elements in the resistivity change matrix represent the resistivity ratio of the resistivity model at adjacent times at each grid cell;

[0012] Determine the first weight coefficient corresponding to each grid cell according to the matrix values of the resistivity change matrix.

[0013] In an optional embodiment, the resistivity change matrix is expressed as:

[0014]

[0015] Among them, represents the resistivity change matrix, represents the initial resistivity model, represents the change in resistivity at the t-th observation time relative to the previous observation time, represents the resistivity ratio of the resistivity model at adjacent times at each grid cell, represents the number of observation data times.

[0016] In an optional embodiment, the step of determining the first weight coefficient corresponding to each grid cell according to the matrix values of the resistivity change matrix includes:

[0017] If the matrix value of the resistivity change matrix is less than a first set value or the matrix value of the resistivity change matrix is greater than a fourth set value, then set the first weight coefficient corresponding to each grid cell to 0;

[0018] If the matrix value of the resistivity change matrix is greater than the first set value and less than the second set value, or the matrix value of the resistivity change matrix is greater than the third set value and less than the fourth set value, then set the first weight coefficient corresponding to each grid cell to a first target value;

[0019] If the matrix value of the resistivity change matrix is greater than the second set value and less than the third set value, set the first weight coefficient corresponding to each grid cell to the second target value;

[0020] Wherein, the first set value is less than the second set value, the second set value is less than the third set value, the third set value is less than the fourth set value, and the first target value is less than the second target value.

[0021] In an alternative embodiment, the step of dividing the research area corresponding to the underground medium into a plurality of grid cells includes:

[0022] Obtain the X-direction data range, Y-direction data range, and Z-direction data range of the observation area;

[0023] Set a target multiple for the X-direction data range, the Y-direction data range, and the Z-direction data range to obtain the data ranges of the research area in the X-direction, Y-direction, and Z-direction; wherein, the target multiple is between 3 times and 5 times.

[0024] Divide into a plurality of grid cells according to the data ranges of the research area in the X-direction, Y-direction, and Z-direction. The plurality of grid cells include non-boundary grid cells and boundary grid cells. The boundary grid cells include boundary surface grid cells, boundary line grid cells, and boundary point grid cells.

[0025] In an alternative embodiment, the second weight coefficient is determined through the following steps:

[0026] For non-boundary grid cells, construct a resistivity change relation formula; determine the second weight coefficient corresponding to each non-boundary grid cell according to the variable result value of the resistivity change relation formula; wherein, the non-boundary grid cell has an upper grid cell, a lower grid cell, a left grid cell, a right grid cell, a front grid cell, and a rear grid cell;

[0027] For boundary grid cells, set the second weight coefficient corresponding to each boundary grid cell to 0.

[0028] In an alternative embodiment, the resistivity change relation formula is expressed as:

[0029]

[0030] Wherein, represents the change amount of the resistivity value of the resistivity model at the non-boundary grid cell at the current moment relative to the resistivity value at the adjacent grid cell, represents the resistivity value of the resistivity model at the non-boundary grid cell at the current moment, represents the resistivity value of the resistivity model at the right grid cell at the current moment, represents the resistivity value of the resistivity model at the front grid cell at the current moment, represents the resistivity value of the resistivity model at the lower grid cell at the current moment, represents the resistivity value of the resistivity model at the left grid cell at the current moment, represents the resistivity value of the resistivity model at the rear grid cell at the current moment, represents the resistivity value of the resistivity model at the upper grid cell at the current moment.

[0031] In an alternative embodiment, determining the second weight coefficient corresponding to each non-boundary grid cell according to the variable result value of the resistivity change relation includes:

[0032] If the variable result value is greater than the first variable threshold, set the second weight coefficient corresponding to each non-boundary grid cell to the third target value;

[0033] If the variable result value is greater than the second variable threshold and not greater than the first variable threshold, set the second weight coefficient corresponding to each non-boundary grid cell to the fourth target value;

[0034] If the variable result value is not greater than the second variable threshold, set the second weight coefficient corresponding to each non-boundary grid cell to the fifth target value;

[0035] Wherein, the first variable threshold is greater than the second variable threshold, the third target value is greater than the fourth target value, and the fourth target value is greater than the fifth target value.

[0036] In an alternative embodiment, the objective function further includes a data term, a model term, and a second weight factor corresponding to the model term. The data term is used to calculate the difference between the apparent resistivity data and the resistivity data corresponding to the inversion result response. The model term is used to calculate the smoothness of the resistivity model. The second weight factor characterizes the proportion of the model term in the objective function.

[0037] In an alternative embodiment, the objective function is represented by the following formula:

[0038]

[0039] Wherein, represents the resistivity model, represents the objective function of time-lapse inversion, represents the data term, represents the model term, represents the time-lapse term, Denote the second weight factor corresponding to the model term in the objective function. Denote the first weight factor corresponding to the time shift term in the objective function. Denote the first weight coefficient. Denote the second weight coefficient. Is the symbol of the Hadamard product, indicating the multiplication of the first weight coefficient and the second weight coefficient corresponding to the same grid cell.

[0040] In a second aspect, an embodiment of the present application further provides a time-lapse inversion device based on time shift term constraints. The time-lapse inversion device includes:

[0041] A data acquisition module, configured to acquire apparent resistivity data of the underground medium.

[0042] A grid division module, configured to divide the research area corresponding to the underground medium into multiple grid cells. The initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model. The spatial range of the research area is larger than the spatial range of the observation area of the underground medium.

[0043] A function construction module, configured to construct an objective function for time-lapse inversion. The objective function includes a time shift term and a first weight factor corresponding to the time shift term. The first weight factor is used to adjust the proportion of the time shift term in the objective function. Wherein, the first weight factor is the product of a first weight coefficient and a second weight coefficient. The first weight coefficient is used to reflect the resistivity change amount of the resistivity model at adjacent times at each grid cell. The second weight coefficient is used to reflect the change amount of the resistivity value of the resistivity model at the current time at each grid cell relative to the resistivity value at the adjacent grid cell.

[0044] A model update module, configured to use the apparent resistivity data to continuously iteratively update the model parameters of the resistivity model to continuously reduce the objective function, so as to obtain an optimal resistivity model.

[0045] In a third aspect, an embodiment of the present application further provides an electronic device, including: a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the electronic device runs, the processor communicates with the memory through the bus. When the machine-readable instructions are executed by the processor, the steps of the time-lapse inversion method as described above are executed.

[0046] In a fourth aspect, an embodiment of the present application further provides a computer-readable storage medium. A computer program is stored on the computer-readable storage medium. When the computer program is run by a processor, the steps of the time-lapse inversion method as described above are executed.

[0047] The embodiment of the present application provides a time-lapse inversion method based on time-lapse term constraint. In this time-lapse inversion method, first, apparent resistivity data of the underground medium is obtained; then, the research area corresponding to the underground medium is divided into multiple grid cells, and the initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model. The spatial range of the research area is larger than the spatial range of the observation area of the underground medium; then, a target function for time-lapse inversion is constructed. The target function includes a time-lapse term and a first weight factor corresponding to the time-lapse term. The first weight factor is used to adjust the proportion of the time-lapse term in the target function. Among them, the first weight factor is the product of a first weight coefficient and a second weight coefficient. The first weight coefficient is used to reflect the resistivity change amount of the resistivity model at each grid cell at adjacent times, and the second weight coefficient is used to reflect the change amount of the resistivity value of the resistivity model at each grid cell at the current time relative to the resistivity value at the adjacent grid cell; finally, using the apparent resistivity data, the model parameters of the resistivity model are continuously iteratively updated to continuously reduce the target function to obtain an optimal resistivity model. The embodiment of the present application selects appropriate weight factors for the time-lapse term, specifically, selects appropriate first weight coefficient and second weight coefficient according to different change situations of the underground area to better reflect the real changes of the underground model, can reduce the error of the inversion result, and thus improve the accuracy of the inversion result.

[0048] To make the above objects, features, and advantages of the present application more obvious and understandable, the following specifically gives preferred embodiments and cooperates with the attached drawings to make a detailed description as follows. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] To more clearly illustrate the technical solutions of the embodiments of the present application, the following will briefly introduce the drawings required to be used in the embodiments. It should be understood that the following drawings only show some embodiments of the present application, so they should not be regarded as limiting the scope. For those of ordinary skill in the art, without creative efforts, other related drawings can also be obtained based on these drawings.

[0050] Figure 1 It is a flowchart of a time-lapse inversion method based on time-lapse term constraint provided by the embodiment of the present application;

[0051] Figure 2 It is a schematic diagram of simulating the change of the underground medium over time in three different situations set by the embodiment of the present application;

[0052] Figure 3 It is a schematic diagram of the division of a grid cell provided by the embodiment of the present application;

[0053] Figure 4 It is a schematic diagram of the structure of a boundary grid cell and a non-boundary grid cell provided by the embodiment of the present application;

[0054] Figure 5 A comparison chart of resistivity inversion results between the traditional method provided by the embodiments of the present application and the method of the present invention;

[0055] Figure 6 A schematic structural diagram of a time-lapse inversion device based on time-lapse term constraint provided by the embodiments of the present application;

[0056] Figure 7 A schematic structural diagram of an electronic device provided by the embodiments of the present application. Specific embodiments

[0057] To make the objectives, technical solutions, and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Usually, the components of the embodiments of the present application described and illustrated in the accompanying drawings here can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present application provided in the accompanying drawings is not intended to limit the scope of the present application to be protected, but only represents the selected embodiments of the present application. Based on the embodiments of the present application, every other embodiment obtained by those skilled in the art without creative efforts belongs to the scope of protection of the present application.

[0058] First, the applicable application scenarios of the present application are introduced. The present application can be applied to the field of geophysical exploration technology. Geophysics is one of the main disciplines of earth science and is a comprehensive discipline for studying the earth and searching for mineral resources inside the earth. According to the different physical fields studied, it can be divided into gravity exploration, magnetic method exploration, electrical method exploration, and seismic exploration, etc. The resistivity method is one of the main branch methods of electrical method exploration. Due to its advantages of simple operation, low cost, high work efficiency, etc., it is widely used in scenarios such as environmental protection, ecological restoration, seawater intrusion, and geological disasters.

[0059] Since the changes in the external environment during the monitoring process affect the monitoring data, resulting in the inversion results of the monitoring data being unable to accurately invert the changes in the underground medium, the time-lapse inversion method based on time-lapse term constraint is often used in the prior art for inversion. Among them, the time-lapse term synchronously inverts the data at different times, and the data at different times constrain each other, which can ensure that the inversion results at different times change continuously with time and has a certain effect of suppressing background noise interference. However, in previous studies, the weight factor of the time-lapse term was usually set as a constant. This method ignores the fact that some areas underground have changed greatly within a certain period of time, rather than changing continuously with time, which easily leads to certain errors in the inversion results and makes the accuracy of the inversion results relatively low.

[0060] Based on this, the embodiments of the present application provide a time-lapse inversion method based on time-lapse term constraints. By selecting appropriate weight factors for the time-lapse terms, the error of the inversion result is reduced, thereby improving the accuracy of the inversion result.

[0061] Please refer to Figure 1 , Figure 1 which is a flowchart of a time-lapse inversion method based on time-lapse term constraints provided by the embodiments of the present application. As Figure 1 shown in

[0062] S101. Obtain the apparent resistivity data of the underground medium;

[0063] S102. Divide the research area corresponding to the underground medium into multiple grid cells. The initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model. The spatial range of the research area is larger than the spatial range of the observation area of the underground medium;

[0064] S103. Construct an objective function for time-lapse inversion. The objective function includes a time-lapse term and a first weight factor corresponding to the time-lapse term. The first weight factor is used to adjust the proportion of the time-lapse term in the objective function. Among them, the first weight factor is the product of a first weight coefficient and a second weight coefficient. The first weight coefficient is used to reflect the resistivity change amount at each grid cell of the resistivity model at adjacent times, and the second weight coefficient is used to reflect the change amount of the resistivity value of the resistivity model at each grid cell at the current time relative to the resistivity value at the adjacent grid cell;

[0065] S104. Use the apparent resistivity data to continuously iterate and update the model parameters of the resistivity model to continuously reduce the objective function, so as to obtain the optimal resistivity model.

[0066] In the above steps S101 to S104, select appropriate weight factors for the time-lapse terms, specifically, select appropriate first weight coefficients and second weight coefficients according to different change situations in the underground area, so as to better reflect the real changes of the underground model, reduce the error of the inversion result, and thus improve the accuracy of the inversion result.

[0067] The following uses specific embodiments to exemplarily illustrate the above steps S101 to S104:

[0068] In step S101, obtain the apparent resistivity data of the underground medium.

[0069] In an alternative embodiment, apparent resistivity data of the subsurface medium can be obtained by arranging survey lines in the subsurface medium; in another alternative embodiment, apparent resistivity data of the subsurface medium can be obtained by means of a flume experiment. Optionally, the apparent resistivity data includes receiver point coordinates, and the spatial range corresponding to the observation area of the subsurface medium is determined by these receiver point coordinates.

[0070] Taking the acquisition of apparent resistivity data of the subsurface medium through a flume experiment as an example for illustration:

[0071] First, the device is arranged: The purpose of time-lapse inversion is to monitor the change process of the subsurface medium. To simulate the change situation of the subsurface medium, copper plates of different specifications are used to simulate the change process of the subsurface target over time. The copper plate is a low-resistivity and high-polarization anomaly with respect to water. Exemplarily, the small copper plate has a length, width, and height of 0.5 m, 0.25 m, and 0.01 m respectively, and the large copper plate has a length, width, and height of 0.6 m, 0.25 m, and 0.03 m respectively. During the data acquisition process, steel electrodes are selected as the power supply electrodes, and calomel electrodes with more stable electrochemical properties are selected as the receiver electrodes.

[0072] Then, data is acquired: The center point of the flume is set as the origin, and the area around the origin is selected as the research area. Figure 2 This is a schematic diagram for simulating the change of the subsurface medium over time in three different situations set in the embodiments of the present application. At time T1, a small copper plate is selected as the anomaly, at time T2, a large copper plate is selected as the anomaly, and at time T3, a combination of small and large copper plates is selected as the anomaly. The depth of the copper plates into the water is 0.15 m. The measuring point positions and observation methods are the same at the three times, and the symmetric four-pole device is used for measurement. The total number of measuring points at a single time is 170, and each measuring point corresponds to 12 different ABs. The spacing information of AB and MN is shown in Table 1; during the acquisition, the power supply period is 8 s, and each data point is iterated at least 5 times.

[0073]

[0074] Table 1. Statistical table of the spacings of AB and MN corresponding to the measuring points in the flume experiment

[0075] In step S102, the research area corresponding to the subsurface medium is divided into multiple grid cells, and the initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model. The spatial range of the research area is larger than the spatial range of the observation area of the subsurface medium.

[0076] In this step, when dividing the research area corresponding to the subsurface medium into multiple grid cells, the research area corresponding to the subsurface medium can be regarded as a cuboid spatial area. As Figure 3 shown, the research area of the subsurface medium can be divided into 64 * 24 * 20 grid cells.

[0077] Among them, the research area is divided into grids mainly to facilitate the mathematical description and numerical calculation of the underground medium. Since the underground medium is a continuous three-dimensional region, it is very difficult to directly analyze and calculate it. Through grid division, the continuous underground medium is discretized into many small grid cells, so that each grid cell can be analyzed and processed separately, thus simplifying the calculation process and helping to improve the accuracy of inversion.

[0078] Here, a homogeneous half-space can be selected as the initial resistivity model. After establishing the initial resistivity model, the initial resistivity value can be assigned to each grid cell according to this model.

[0079] That is to say, under the condition of taking the homogeneous half-space as the initial resistivity model, when the grid cells are divided, the initial resistivity value of each grid cell is set to the same value. With the acquisition of apparent resistivity data, this initial homogeneous half-space model and the resistivity values of each grid cell can be adjusted and corrected to more accurately reflect the actual underground resistivity distribution.

[0080] In an optional embodiment, step S102 specifically includes:

[0081] Step S1021: Obtain the data ranges in the X direction, Y direction, and Z direction of the observation area.

[0082] Here, the data ranges in the X direction, Y direction, and Z direction of the observation area can be determined according to the coordinate ranges in the X direction, Y direction, and Z direction of the receiving points. It can also be understood that the coordinate ranges in the X direction, Y direction, and Z direction of the receiving points are the data ranges in the X direction, Y direction, and Z direction of the observation area.

[0083] Step S1022: Set target multiples for the data ranges in the X direction, Y direction, and Z direction to obtain the data ranges of the research area in the X direction, Y direction, and Z direction; where the target multiple is between 3 times and 5 times.

[0084] In the embodiment of the present application, if the entire spatial range of the research area is set too small, boundary effects will occur during the inversion process and the inversion result will be poor. Therefore, the data range of the research area is generally set to 3 - 5 times the data range of the observation area, which can obtain a more accurate inversion result, thereby reducing the boundary effects during the inversion process and improving the accuracy of the inversion result.

[0085] Step S1023: Divide the data range of the research area in the X, Y, and Z directions into multiple grid cells, where the multiple grid cells include non-boundary grid cells and boundary grid cells, and the boundary grid cells include surface boundary grid cells, boundary line grid cells, and boundary point grid cells.

[0086] Exemplarily, as Figure 4 shown, (a), (b), (c), and (d) respectively show non-boundary grid cells, surface boundary grid cells, boundary line grid cells, and boundary point grid cells. Taking the non-boundary grid cell c in (a) as an example, the subsequent i, j, and k represent the positions where the grid cell c is located, that is, i, j, and k respectively represent the grid order in the X, Y, and Z directions. Since c is a non-boundary grid cell, it has a right grid cell r, a front grid cell p, a bottom grid cell b, a left grid cell l, a back grid cell n, and an upper grid cell u.

[0087] Specifically, since the underground medium is divided into three-dimensional grid cuboids, for (b), (c), and (d), they respectively correspond to the surface boundary grid cells (cells on the cuboid surface), boundary line grid cells (cells on the cuboid edges), and boundary point grid cells (cells at the cuboid vertices) of the grid cuboid.

[0088] In step S103, a time-lapse inversion objective function is constructed. The objective function includes a time-lapse term and a first weight factor corresponding to the time-lapse term. The first weight factor is used to adjust the proportion of the time-lapse term in the objective function. Among them, the first weight factor is the product of a first weight coefficient and a second weight coefficient. The first weight coefficient is used to reflect the resistivity change amount of the resistivity model at each grid cell at adjacent times, and the second weight coefficient is used to reflect the change amount of the resistivity value of the resistivity model at each grid cell at the current time relative to the resistivity value at the adjacent grid cell.

[0089] Here, multiplying the first weight coefficient corresponding to each grid cell by the second weight coefficient can obtain the first weight factor corresponding to the time-lapse term. Specifically, the resistivity change amount of the resistivity model at each grid cell at adjacent times can be reflected by the first weight coefficient, and the change amount of the resistivity value of the resistivity model at each grid cell at the current time relative to the resistivity value at the adjacent grid cell can be reflected by the second weight coefficient. Then, the product of the first weight coefficient and the second weight coefficient is used as the first weight factor of the time-lapse term to adjust the proportion of the time-lapse term in the objective function, which can better handle the changes in the underground gradual region and mutation region to reduce the error of time-lapse inversion.

[0090] In an optional embodiment, the first weight coefficient is determined through the following steps:

[0091] Construct a resistivity change matrix, where the matrix elements in the resistivity change matrix represent the resistivity ratios of the resistivity models at adjacent times at each grid cell;

[0092] Determine the first weight coefficient corresponding to each grid cell according to the matrix values of the resistivity change matrix.

[0093] Here, according to the magnitudes of the matrix values of the resistivity change matrix, different values can be selected for the first weight coefficient of each grid cell. When the matrix value of the resistivity change matrix is too large or too small, a smaller first weight coefficient is selected; otherwise, a larger first weight coefficient is selected. The effect of such a value selection is to allow sudden changes in the physical properties of the underground model, prevent the inversion results at different times from being too close, so as to reflect the real changes of the underground model, and thus obtain a better inversion result.

[0094] Specifically, the resistivity change matrix is expressed as:

[0095]

[0096] Among them, represents the resistivity change matrix, represents the initial resistivity model, represents the resistivity change amount at the t-th observation time relative to the previous observation time, represents the resistivity ratio of the resistivity models at adjacent times at each grid cell, represents the number of observation data times.

[0097] Furthermore, the steps of determining the first weight coefficient corresponding to each grid cell according to the matrix values of the resistivity change matrix include:

[0098] If the matrix value of the resistivity change matrix is less than the first set value or the matrix value of the resistivity change matrix is greater than the fourth set value, set the first weight coefficient corresponding to each grid cell to 0;

[0099] If the matrix value of the resistivity change matrix is greater than the first set value and less than the second set value, or the matrix value of the resistivity change matrix is greater than the third set value and less than the fourth set value, set the first weight coefficient corresponding to each grid cell to the first target value;

[0100] If the matrix value of the resistivity change matrix is greater than the second set value and less than the third set value, set the first weight coefficient corresponding to each grid cell to the second target value;

[0101] Among them, the first set value is less than the second set value, the second set value is less than the third set value, the third set value is less than the fourth set value, and the first target value is less than the second target value.

[0102] Exemplarily, the first set value can be set to 0.8, the second set value can be set to 0.85, the third set value can be set to 1.15, the fourth set value can be set to 1.2, the first target value can be set to 0.01, and the second target value can be set to 1. Specifically, as shown in Table 2:

[0103]

[0104] Table 2. Corresponding relationship between matrix values and the first weight coefficient

[0105] In an alternative embodiment, the second weight coefficient is determined through the following steps:

[0106] For non-boundary grid cells, construct a resistivity change relationship; based on the variable result values of the resistivity change relationship, determine the second weight coefficient corresponding to each non-boundary grid cell; wherein, the non-boundary grid cells have upper grid cells, lower grid cells, left grid cells, right grid cells, front grid cells, and rear grid cells;

[0107] For boundary grid cells, set the second weight coefficient corresponding to each boundary grid cell to 0.

[0108] In an actual scenario, since the grid setting range of the research area is larger than the grid setting range of the actual observation area, for boundary grid cells, the model change rate thereof can be not considered, and thus the second weight coefficients of boundary grid cells are all 0.

[0109] Here, each non-boundary grid cell needs to obtain a variable result value according to the resistivity change relationship, and the variable result values of each non-boundary grid cell are different.

[0110] In an alternative embodiment, the resistivity change relationship is expressed as:

[0111]

[0112] Wherein, represents the change amount of the resistivity value of the resistivity model at the non-boundary grid cell at the current moment relative to the resistivity value at the adjacent grid cell, represents the resistivity value of the resistivity model at the non-boundary grid cell at the current moment, represents the resistivity value of the resistivity model at the right grid cell at the current moment, represents the resistivity value of the resistivity model at the front grid cell at the current moment, represents the resistivity value of the resistivity model at the lower grid cell at the current moment, represents the resistivity value of the resistivity model at the left grid cell at the current moment, represents the resistivity value of the resistivity model at the rear grid cell at the current moment, represents the resistivity value of the resistivity model at the upper grid cell at the current moment.

[0113] In an alternative embodiment, the step of determining the second weight coefficient corresponding to each non-boundary grid cell according to the variable result value of the resistivity change relation includes:

[0114] If the variable result value is greater than the first variable threshold, set the second weight coefficient corresponding to each non-boundary grid cell to the third target value;

[0115] If the variable result value is greater than the second variable threshold and not greater than the first variable threshold, set the second weight coefficient corresponding to each non-boundary grid cell to the fourth target value;

[0116] If the variable result value is not greater than the second variable threshold, set the second weight coefficient corresponding to each non-boundary grid cell to the fifth target value;

[0117] wherein, the first variable threshold is greater than the second variable threshold, the third target value is greater than the fourth target value, and the fourth target value is greater than the fifth target value.

[0118] Here, it is possible to judge the resistivity change situation of the resistivity model at each grid cell at the current moment according to If the value of is relatively large, it represents that the resistivity model has changed at this position at the current moment, then judge that this grid cell is a gradual change region, and the second weight coefficient can be set to 1. If the value of is medium, then the second weight coefficient can be set to 0.1; if the value of is too small, it represents that the resistivity model has no obvious change at this position at the current moment, then judge that this grid cell is a mutation region (regardless of whether there is a change at the next moment, the change situation at this grid cell does not need to be considered), and the second weight coefficient can be set to 0.01.

[0119] Exemplarily, the first variable threshold can be set to 0.5, the second variable threshold can be set to 0.2, the third target value can be set to 1, the fourth target value can be set to 0.1, and the fifth target value can be set to 0.01. Specifically, as shown in Table 3:

[0120]

[0121] Table 3. Corresponding relationship between variable result value and second weight coefficient

[0122] In an alternative embodiment, the objective function further includes a data term, a model term, and a second weight factor corresponding to the model term. The data term is used to calculate the difference between the apparent resistivity data and the resistivity data corresponding to the inversion result response. The model term is used to calculate the smoothness of the resistivity model. The second weight factor represents the proportion of the model term in the objective function.

[0123] Here, the data term, the model term, and the second weight factor cooperate with each other in time-lapse inversion and jointly act on the objective function to guide the optimization and adjustment of the resistivity model, ultimately obtaining a more accurate and reasonable subsurface model.

[0124] Specifically, the objective function can be expressed by the following formula:

[0125]

[0126] Where, represents the resistivity model, represents the objective function of time-lapse inversion, represents the data term, represents the model term, represents the time-lapse term, represents the second weight factor corresponding to the model term in the objective function, represents the first weight factor corresponding to the time-lapse term in the objective function, represents the first weight coefficient, represents the second weight coefficient, is the Hadamard product symbol, indicating the multiplication of the first weight coefficient and the second weight coefficient corresponding to the same grid cell.

[0127] Here, the inversion is achieved by continuously minimizing the objective function, will change continuously, and finally the resistivity model that best conforms to the subsurface situation is obtained.

[0128] In step S104, using the apparent resistivity data, the model parameters of the resistivity model are continuously iteratively updated to continuously reduce the objective function, so as to obtain the optimal resistivity model.

[0129] In this step, after multiple iterative updates, the objective function is minimized. When certain convergence conditions are met (such as the objective function value no longer decreases significantly, the parameter change amount is less than a certain threshold, etc.), the final resistivity model is obtained. The resistivity model describes the distribution of the resistivity of the subsurface medium. These models can visually display the electrical properties at different positions underground and provide important reference bases for geological interpretation, resource exploration (such as searching for mineral resources, groundwater resources, etc.), engineering construction (such as evaluating the foundation stability, etc.).

[0130] For example, during the experiment, a small copper plate was selected as the anomaly at time T1, that is, a small copper plate was inserted on the left side of the water tank. A large copper plate was selected as the anomaly at time T2, that is, a large copper plate was inserted on the left side of the water tank. A combination of large and small copper plates was selected as the anomaly at time T3, that is, a large copper plate and a small copper plate were inserted on the left and right sides of the water tank respectively. Under the above experimental conditions, the Figure 5 inversion results shown can be obtained. Specifically, Figure 5 FIG. Figure 5 is a comparison chart of the resistivity inversion results of the traditional method and the method of the present invention provided by the embodiments of the present application, showing the resistivity inversion result profile at a depth of Z = 0.15 m. Among them, the white marked circles and black marked circles are used to mark the anomaly positions in the inversion results. It should be noted that since the black marked circles cannot be shown on the black background, the white marked circles are used instead. It can be seen that in the results of the traditional method (left), anomalies appear at both the left and right positions at all three times, and the anomaly position at the third time has a large difference from the true position of the copper plate; the anomaly positions in the method of the present invention (right) are all close to the true positions of the copper plates. That is to say, when using the traditional method for time-lapse inversion, anomalies also appear on the right side at the first two times but there is actually no copper plate on the right side. At the same time, the anomaly position at the third time has a large difference from the true position of the copper plate. The above anomalies are all caused by improper setting of the weight factor during the time-lapse inversion process. It can be seen that the method of the present invention can avoid this problem and has obvious advantages.

[0131] The time-lapse inversion method based on time-lapse term constraint provided by the embodiments of the present application selects appropriate weight factors for the time-lapse terms, specifically, selects appropriate first and second weight coefficients according to different change situations of the underground area to better reflect the true changes of the underground model, can reduce the error of the inversion results, and thus improve the accuracy of the inversion results.

[0132] Based on the same inventive concept, the embodiments of the present application also provide a time-lapse inversion device based on time-lapse term constraint corresponding to the time-lapse inversion method based on time-lapse term constraint. Since the principle of solving problems by the device in the embodiments of the present application is similar to the time-lapse inversion method based on time-lapse term constraint described above in the embodiments of the present application, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be elaborated.

[0133] Please refer to Figure 6 , Figure 6 which is a schematic structural diagram of a time-lapse inversion device based on time-lapse term constraint provided by the embodiments of the present application. As shown in Figure 6 , the time-lapse inversion device 700 includes:

[0134] A data acquisition module 701, configured to acquire apparent resistivity data of underground media;

[0135] The mesh generation module 702 is configured to divide the research area corresponding to the underground medium into a plurality of grid cells. The initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model. The spatial range of the research area is larger than the spatial range of the observation area of the underground medium;

[0136] The function construction module 703 is configured to construct an objective function for time-lapse inversion. The objective function includes a time-lapse term and a first weight factor corresponding to the time-lapse term. The first weight factor is used to adjust the proportion of the time-lapse term in the objective function. Among them, the first weight factor is the product of a first weight coefficient and a second weight coefficient. The first weight coefficient is used to reflect the resistivity change amount of the resistivity model at adjacent times at each grid cell, and the second weight coefficient is used to reflect the change amount of the resistivity value of the resistivity model at the current time at each grid cell relative to the resistivity value at an adjacent grid cell;

[0137] The model update module 704 is configured to continuously iteratively update the model parameters of the resistivity model by using the apparent resistivity data to continuously reduce the objective function, so as to obtain an optimal resistivity model.

[0138] Please refer to Figure 7 , Figure 7 which is a schematic structural diagram of an electronic device provided by an embodiment of the present application. As shown in Figure 7 , the electronic device 800 includes a processor 801, a memory 802, and a bus 803.

[0139] The memory 802 stores machine-readable instructions executable by the processor 801. When the electronic device 800 runs, the processor 801 communicates with the memory 802 through the bus 803. When the machine-readable instructions are executed by the processor 801, the steps of the time-lapse inversion method based on time-lapse term constraint in the method embodiment as shown above can be executed. The specific implementation manner can refer to the method embodiment and will not be elaborated here. Figure 1 shown.

[0140] An embodiment of the present application further provides a computer-readable storage medium. A computer program is stored on the computer-readable storage medium. When the computer program is run by a processor, the steps of the time-lapse inversion method based on time-lapse term constraint in the method embodiment as shown above can be executed. The specific implementation manner can refer to the method embodiment and will not be elaborated here. Figure 1 shown.

[0141] Those skilled in the art can clearly understand that for the convenience and conciseness of description, the specific working processes of the systems, devices, and units described above can refer to the corresponding processes in the foregoing method embodiments and will not be elaborated here.

[0142] In several embodiments provided in the present application, it should be understood that the disclosed systems, devices, and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of the units is only a logical function division. In actual implementation, there may be other division methods. For another example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the displayed or discussed coupling or direct coupling or communication connection between each other can be through some communication interfaces. The indirect coupling or communication connection of the devices or units can be in electrical, mechanical, or other forms.

[0143] The units described as separate components may or may not be physically separated. The components displayed as units may or may not be physical units, that is, they can be located in one place, or they can be distributed to multiple network units. Some or all of the units can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0144] In addition, in each embodiment of the present application, the functional units can be integrated in a processing unit, or each unit can exist physically alone, or two or more units can be integrated in one unit.

[0145] If the function is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a non-volatile computer-readable storage medium executable by a processor. Based on such an understanding, the technical solution of the present application, in essence, or the part that contributes to the prior art, or this part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in each embodiment of the present application. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical discs that can store program codes.

[0146] Finally, it should be noted that the above-described embodiments are only specific implementation manners of the present application, used to illustrate the technical solutions of the present application, rather than limiting it. The protection scope of the present application is not limited thereto. Although the present application has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that any person skilled in the art within the technical scope disclosed by the present application can still modify the technical solutions recorded in the foregoing embodiments, or can easily think of changes, or perform equivalent replacements on some of the technical features; and these modifications, changes or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present application, and should all be covered within the protection scope of the present application. Therefore, the protection scope of the present application shall be subject to the protection scope of the claims.

Claims

1. A time-shift inversion method based on time-shift term constraints, characterized in that: The time-shift inversion method comprises: Obtain apparent resistivity data of underground media; Dividing the study area corresponding to the underground medium into a plurality of grid cells, the initial resistivity value of each grid cell is determined according to a pre-established initial resistivity model, and the spatial range of the study area is larger than the spatial range of the observation area of ​​the underground medium; Constructing an objective function of time-shift inversion, the objective function comprising a time-shift term and a first weight factor corresponding to the time-shift term, the first weight factor being used to adjust the proportion of the time-shift term in the objective function, wherein the first weight factor is the product of a first weight coefficient and a second weight coefficient, the first weight coefficient being used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient being used to reflect the resistivity value of the resistivity model at each grid unit at the current moment relative to the resistivity value at the adjacent grid unit; The apparent resistivity data is used to continuously iteratively update the model parameters of the resistivity model so as to continuously reduce the objective function, thereby obtaining an optimal resistivity model.

2. The time-lapse inversion method according to claim 1, characterized in that: The first weight coefficient is determined by the following steps: Constructing a resistivity variation matrix, wherein the matrix elements in the resistivity variation matrix represent the resistivity ratio of the resistivity model at each grid unit at adjacent moments; A first weight coefficient corresponding to each grid unit is determined according to the matrix value of the resistivity variation matrix.

3. The time-lapse inversion method according to claim 2, characterized in that: The resistivity variation matrix is ​​expressed as: in, represents the resistivity variation matrix, represents the initial resistivity model, represents the resistivity change at the tth observation time relative to the previous observation time, represents the resistivity ratio of the resistivity model at each grid unit at adjacent moments, Indicates the number of observations.

4. The time-lapse inversion method according to claim 2, characterized in that: Determining a first weight coefficient corresponding to each grid unit according to the matrix value of the resistivity variation matrix includes: If the matrix value of the resistivity change matrix is ​​less than the first set value or the matrix value of the resistivity change matrix is ​​greater than the fourth set value, setting the first weight coefficient corresponding to each grid unit to 0; If the matrix value of the resistivity change matrix is ​​greater than the first set value and less than the second set value, or the matrix value of the resistivity change matrix is ​​greater than the third set value and less than the fourth set value, then the first weight coefficient corresponding to each grid unit is set to the first target value; If the matrix value of the resistivity change matrix is ​​greater than the second set value and less than the third set value, the first weight coefficient corresponding to each grid unit is set to the second target value; Among them, the first set value is smaller than the second set value, the second set value is smaller than the third set value, the third set value is smaller than the fourth set value, and the first target value is smaller than the second target value.

5. The time-lapse inversion method according to claim 1, characterized in that: The step of dividing the research area corresponding to the underground medium into a plurality of grid units includes: Obtaining the data range in the X direction, the data range in the Y direction, and the data range in the Z direction of the observation area; Setting a target multiple for the X-direction data range, the Y-direction data range, and the Z-direction data range to obtain the data range of the study area in the X-direction, the Y-direction, and the Z-direction; wherein the target multiple is between 3 and 5 times; According to the data range of the study area in the X direction, the Y direction and the Z direction, the study area is divided into a plurality of grid cells, the plurality of grid cells include non-boundary grid cells and boundary grid cells, and the boundary grid cells include boundary surface grid cells, boundary line grid cells and boundary point grid cells.

6. The time-lapse inversion method according to claim 5, characterized in that: The second weight coefficient is determined by the following steps: For non-boundary grid cells, a resistivity change relationship is constructed; according to the variable result value of the resistivity change relationship, a second weight coefficient corresponding to each non-boundary grid cell is determined; wherein the non-boundary grid cell includes an upper grid cell, a lower grid cell, a left grid cell, a right grid cell, a front grid cell and a rear grid cell; For the boundary grid cells, the second weight coefficient corresponding to each boundary grid cell is set to 0.

7. The time-lapse inversion method according to claim 6, characterized in that: The resistivity change relationship is expressed as: in, It represents the change of the resistivity value of the resistivity model at the non-boundary grid unit relative to the resistivity value at the adjacent grid unit at the current moment. represents the resistivity value of the resistivity model at the non-boundary grid unit at the current moment, Indicates the resistivity value of the resistivity model at the right grid unit at the current moment, represents the resistivity value of the resistivity model at the front grid unit at the current moment, represents the resistivity value of the resistivity model at the lower grid unit at the current moment, Indicates the resistivity value of the resistivity model at the left grid cell at the current moment, represents the resistivity value of the resistivity model at the rear grid unit at the current moment, Represents the resistivity value of the resistivity model at the upper grid cell at the current moment.

8. The time-lapse inversion method according to claim 7, characterized in that: Determining the second weight coefficient corresponding to each non-boundary grid unit according to the variable result value of the resistivity change relationship includes: If the variable result value is greater than the first variable threshold, the second weight coefficient corresponding to each non-boundary grid unit is set to a third target value; If the variable result value is greater than the second variable threshold value and not greater than the first variable threshold value, setting the second weight coefficient corresponding to each non-boundary grid unit to a fourth target value; If the variable result value is not greater than the second variable threshold, setting the second weight coefficient corresponding to each non-boundary grid unit to the fifth target value; Among them, the first variable threshold is greater than the second variable threshold, the third target value is greater than the fourth target value, and the fourth target value is greater than the fifth target value.

9. The time-lapse inversion method according to claim 1, characterized in that: The objective function also includes a data item, a model item and a second weight factor corresponding to the model item, the data item is used to calculate the difference between the apparent resistivity data and the resistivity data of the inversion result response, the model item is used to calculate the smoothness of the resistivity model, and the second weight factor represents the proportion of the model item in the objective function.

10. The time-lapse inversion method according to claim 9, characterized in that: The objective function is expressed by the following formula: in, represents the resistivity model, represents the objective function of time-lapse inversion, Represents a data item, represents the model term, represents the time-shift term, represents the second weight factor corresponding to the model term in the objective function, represents the first weight factor corresponding to the time shift term in the objective function, represents the first weight coefficient, represents the second weight coefficient, is the Hadamard product symbol, indicating that the first weight coefficient and the second weight coefficient corresponding to the same grid unit are multiplied.

Images