Time shifting inversion method based on time shifting item constraint

Abstract

The invention provides a time shifting inversion method based on time shifting item constraint. The method comprises the following steps: acquiring apparent resistivity data; a research area of an underground medium is divided into a plurality of grid units, an initial resistivity value of each grid unit is determined according to an initial resistivity model, an objective function is constructed, the objective function comprises a time shifting item and a first weight factor, the first weight factor is a product of a first weight coefficient and a second weight coefficient, and the second weight coefficient is a product of a second weight coefficient and a third weight coefficient; the first weight coefficient is used for reflecting the resistivity variation of the resistivity model at each grid unit at the adjacent moment, and the second weight coefficient is used for reflecting the variation of the resistivity value of the resistivity model at each grid unit at the current moment relative to the resistivity value of the resistivity model at the adjacent grid unit; and continuously iterating the model parameters of the resistivity model by using the apparent resistivity data so as to continuously reduce the target function, thereby obtaining the optimal resistivity model. According to the method, the appropriate weight factor is selected for the time shifting item, so that the inversion result error is reduced, and the precision of the inversion result is improved.

Description

Technical Field

[0001] The present application relates to the field of geophysical exploration technology, and in particular to a time-shift inversion method based on time-shift term constraints. Background Art

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

[0003] Since changes in the external environment during the monitoring process affect the monitoring data, the inversion results of the monitoring data cannot accurately invert the changes in the underground medium. In the prior art, a time-shift inversion method based on time-shift term constraints is often used for inversion. The time-shift term inverts data at different times synchronously, and the data at different times constrain each other, which can ensure that the inversion results at different times change continuously over time, and has a certain effect of suppressing background noise interference. However, in previous studies, the weight factor of the time-shift term is usually set to a constant. This method ignores the fact that some underground areas have undergone large changes within a certain period of time, rather than continuously changing over time, which easily leads to certain errors in the inversion results, resulting in low accuracy of the inversion results. Summary of the invention

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

[0005] In a first aspect, an embodiment of the present application provides a time-shift inversion method based on a time-shift term constraint, the time-shift inversion method comprising: 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.

[0006] In an optional embodiment, 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.

[0007] In an optional embodiment, the resistivity variation matrix is ​​expressed as:

[0008] 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.

[0009] In an optional embodiment, determining the 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.

[0010] In an optional embodiment, dividing the study 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.

[0011] In an optional embodiment, 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.

[0012] In an optional embodiment, the resistivity change relationship is expressed as:

[0013] 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.

[0014] In an optional embodiment, 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.

[0015] In an optional embodiment, 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 characterizes the proportion of the model item in the objective function.

[0016] In an optional embodiment, the objective function is expressed by the following formula:

[0017] 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.

[0018] In a second aspect, an embodiment of the present application further provides a time-shift inversion device based on time-shift term constraints, the time-shift inversion device comprising: A data acquisition module, used to acquire apparent resistivity data of underground media; A grid division module, used for dividing the study area corresponding to the underground medium into a plurality of grid units, wherein the initial resistivity value of each grid unit 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; A function construction module, used to construct an objective function of time-shift inversion, wherein the objective function includes a time-shift term and a first weight factor corresponding to the time-shift term, wherein 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, wherein the first weight coefficient is used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient is 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 model updating module is used to utilize the apparent resistivity data to continuously and iteratively update the model parameters of the resistivity model so as to continuously reduce the objective function and obtain an optimal resistivity model.

[0019] In a third aspect, an embodiment of the present application further provides an electronic device, comprising: a processor, a memory and a bus, wherein the memory stores machine-readable instructions executable by the processor, and when the electronic device is running, the processor and the memory communicate through the bus, and when the machine-readable instructions are executed by the processor, the steps of the time-shift inversion method as described above are performed.

[0020] In a fourth aspect, an embodiment of the present application further provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the time-shift inversion method as described above are executed.

[0021] The embodiment of the present application provides a time-shift inversion method based on time-shift term constraints. In the time-shift inversion method, first, the apparent resistivity data of the underground medium is obtained; then, the study area corresponding to the underground medium is divided into multiple grid units, the initial resistivity value of each grid unit 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; then, an objective function of the time-shift inversion is constructed, 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 the first weight coefficient and the second weight coefficient, the first weight coefficient is used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient is 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; finally, using the apparent resistivity data, the model parameters of the resistivity model are continuously iteratively updated to continuously reduce the objective function to obtain the optimal resistivity model. The embodiment of the present application selects a suitable weight factor for the time-shift term, specifically, selects a suitable first weight coefficient and a second weight coefficient according to different changes 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.

[0022] In order to make the above-mentioned objects, features and advantages of the present application more obvious and easy to understand, preferred embodiments are specifically cited below and described in detail with reference to the attached drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings required for use in the embodiments will be briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present application and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without paying creative work.

[0024] Figure 1 A flowchart of a time-shift inversion method based on time-shift term constraints provided in an embodiment of the present application; Figure 2 Schematic diagram of three different situations for simulating the change of underground media over time set in the embodiment of the present application; Figure 3 A schematic diagram of the division of a grid unit provided in an embodiment of the present application; Figure 4 A schematic diagram of the structure of a boundary grid unit and a non-boundary grid unit provided in an embodiment of the present application; Figure 5 A comparison diagram of resistivity inversion results of the conventional method and the method of the present invention provided in the embodiments of the present application; Figure 6 A schematic structural diagram of a time-shift inversion device based on time-shift term constraints provided in an embodiment of the present application; Figure 7 A schematic diagram of the structure of an electronic device provided in an embodiment of the present application. DETAILED DESCRIPTION

[0025] To make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be clearly and completely described below in conjunction with the 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 of the embodiments. The components of the embodiments of the present application usually described and shown in the 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 drawings is not intended to limit the scope of the application claimed for protection, but merely represents the selected embodiments of the present application. Based on the embodiments of the present application, each other embodiment obtained by those skilled in the art without making creative work belongs to the scope of protection of the present application.

[0026] First, the application scenarios to which this application can be applied are introduced. This application can be applied to the field of geophysical exploration technology. Geophysics is one of the main disciplines of earth science and a comprehensive discipline that studies the earth and searches for mineral resources inside the earth. Depending on the physical field studied, it can be divided into gravity exploration, magnetic exploration, electrical exploration, and seismic exploration. Resistivity method is one of the main branch methods of electrical exploration. Due to its advantages such as simple operation, low cost and high work efficiency, it is widely used in scenarios such as environmental protection, ecological restoration, seawater intrusion, and geological disasters.

[0027] Since changes in the external environment during the monitoring process affect the monitoring data, the inversion results of the monitoring data cannot accurately invert the changes in the underground medium. In the prior art, a time-shift inversion method based on time-shift term constraints is often used for inversion. The time-shift term inverts data at different times synchronously, and the data at different times constrain each other, which can ensure that the inversion results at different times change continuously over time, and has a certain effect of suppressing background noise interference. However, in previous studies, the weight factor of the time-shift term is usually set to a constant. This method ignores the fact that some underground areas have undergone large changes within a certain period of time, rather than continuously changing over time, which easily leads to certain errors in the inversion results, resulting in low accuracy of the inversion results.

[0028] Based on this, an embodiment of the present application provides a time-shift inversion method based on time-shift term constraints, which reduces the error of the inversion result by selecting a suitable weight factor for the time-shift term, thereby improving the accuracy of the inversion result.

[0029] See also Figure 1 , Figure 1 The flowchart of a time-shift inversion method based on time-shift term constraints provided in an embodiment of the present application is as follows. Figure 1 As shown in , the time-shift inversion method provided by the embodiment of the present application includes: S101, obtaining apparent resistivity data of underground media; S102, dividing the study area corresponding to the underground medium into a plurality of grid cells, wherein 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; S103, constructing an objective function of time-shift inversion, wherein the objective function includes a time-shift term and a first weight factor corresponding to the time-shift term, wherein 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 the first weight coefficient and the second weight coefficient, wherein the first weight coefficient is used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient is 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; S104, using the apparent resistivity data, continuously iteratively updating the model parameters of the resistivity model so that the objective function is continuously reduced, so as to obtain the optimal resistivity model.

[0030] In the above steps S101 to S104, appropriate weight factors are selected for the time-shift terms, specifically, appropriate first weight coefficients and second weight coefficients are selected according to different changes in the underground area, so as to better reflect the real changes of the underground model, reduce the error of the inversion results, and thus improve the accuracy of the inversion results.

[0031] The above steps S101 to S104 are exemplarily described below through a specific embodiment: In step S101, the apparent resistivity data of the underground medium is obtained.

[0032] In an optional embodiment, the apparent resistivity data of the underground medium can be obtained by arranging a survey line in the underground medium; in another optional embodiment, the apparent resistivity data of the underground medium can be obtained by a water tank experiment method. Optionally, the apparent resistivity data includes the coordinates of the receiving point, and the spatial range corresponding to the observation area of ​​the underground medium is determined by the coordinates of the receiving point.

[0033] Take the water tank experiment to obtain the apparent resistivity data of underground media as an example to illustrate: First, the device is arranged: the purpose of time-lapse inversion is to monitor the change process of underground media. In order to simulate the change of underground media, copper plates of different specifications are used to simulate the change process of underground target bodies over time. Copper plates are low-resistance and high-polarization abnormal bodies relative to water. For example, the length, width and height of small copper plates are 0.5m, 0.25m and 0.01m respectively, and the length, width and height of large copper plates are 0.6m, 0.25m and 0.03m respectively. During data acquisition, steel electrodes are selected as power supply electrodes, and calomel electrodes with more stable electrochemical properties are selected as receiving electrodes.

[0034] Then, data collection is performed: the center point of the water tank is set as the origin, and the area around the origin is selected as the research area. Figure 2 The schematic diagram of the three different situations set for the embodiment of the present application to simulate the change of underground media over time, a small copper plate is selected as the anomaly at time T1, a large copper plate is selected as the anomaly at time T2, and a combination of large and small copper plates is selected as the anomaly at time T3, and the depth of the copper plate into the water is 0.15m. The measurement point positions and observation methods at the three moments are the same, and a symmetrical quadrupole device is used for measurement. The total number of measurement points at a single moment is 170, and each measurement point corresponds to 12 different ABs. The spacing information of AB and MN is shown in Table 1; the power supply cycle during the collection is 8s, and each data point is iterated at least 5 times.

[0035]

[0036] Table 1. Statistics of AB and MN spacing corresponding to the water tank test points In step S102, the study area corresponding to the underground medium is divided into multiple grid units, the initial resistivity value of each grid unit 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.

[0037] In this step, the study area corresponding to the underground medium is divided into multiple grid cells. The study area corresponding to the underground medium can be regarded as a rectangular space area, such as Figure 3 As shown, the study area of ​​the underground medium can be divided into 64*24*20 grid cells.

[0038] The main purpose of dividing the study area into grids is to facilitate the mathematical description and numerical calculation of the underground medium. Since the underground medium is a continuous three-dimensional area, it is very difficult to analyze and calculate it directly. Through grid division, the continuous underground medium is discretized into many small grid units, so that each grid unit can be analyzed and processed separately, which simplifies the calculation process and helps to improve the accuracy of the inversion.

[0039] Here, the uniform half space can be selected as the initial resistivity model. After the initial resistivity model is established, the initial resistivity value can be assigned to each grid unit according to the model.

[0040] That is to say, under the condition of taking the uniform half-space as the initial resistivity model, after the grid units are divided, the initial resistivity value of each grid unit is set to the same value. As the apparent resistivity data is acquired, this initial uniform half-space model and the resistivity value of each grid unit can be adjusted and corrected to more accurately reflect the actual underground resistivity distribution.

[0041] In an optional embodiment, step S102 specifically includes: Step S1021, 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.

[0042] Here, the X-direction data range, Y-direction data range, and Z-direction data range of the observation area can be determined according to the X-direction data range, Y-direction data range, and Z-direction data range of the receiving point coordinates. It can also be understood that the X-direction data range, Y-direction data range, and Z-direction data range of the receiving point coordinates are the X-direction data range, Y-direction data range, and Z-direction data range of the observation area.

[0043] Step S1022, set target multiples for the data range in the X direction, the data range in the Y direction, and the data range in the Z direction 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 times and 5 times.

[0044] In the embodiment of the present application, if the entire spatial range of the study area is set too small, a boundary effect will occur during the inversion process and the inversion result will be poor. Therefore, the data range of the study area is generally set to 3 to 5 times the data range of the observation area. More accurate inversion results can be obtained, thereby reducing the boundary effect in the inversion process and improving the accuracy of the inversion results.

[0045] Step S1023: Divide the study area into a plurality of grid cells according to the data range in the X direction, the Y direction and the 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.

[0046] For example, Figure 4As shown, (a), (b), (c), and (d) respectively show non-boundary grid cells, boundary surface grid cells, boundary line grid cells, and boundary point grid cells. Taking the non-boundary grid cell c in (a) as an example, the following i, j, and k represent the position of the grid cell c, that is, i, j, and k represent the grid order in the three directions of X, Y, and Z, respectively. 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 rear grid cell n, and an upper grid cell u.

[0047] Specifically, The underground medium is divided into three-dimensional grid cuboids. For (b), (c), and (d), the corresponding ones are the boundary surface grid cells (cells on the surface of the cuboid), the boundary line grid cells (cells on the edge of the cuboid), and the boundary point grid cells (cells on the vertices of the cuboid).

[0048] In step S103, an objective function of time-shift inversion is constructed, wherein the objective function includes a time-shift term and a first weight factor corresponding to the time-shift term, wherein 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 the first weight coefficient and the second weight coefficient, wherein the first weight coefficient is used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient is 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.

[0049] Here, the first weight factor corresponding to the time-shift term can be obtained by multiplying the first weight coefficient corresponding to each grid unit with the second weight coefficient. Specifically, the first weight coefficient can be used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient can be used to reflect the change of 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. Then, the product of the first weight coefficient and the second weight coefficient is used as the first weight factor of the time-shift term to adjust the proportion of the time-shift term in the objective function, which can better cope with the changes in the underground gradual change area and the sudden change area to reduce the error of time-shift inversion.

[0050] In an optional embodiment, the first weight coefficient is determined by the following steps: A resistivity variation matrix is ​​constructed, wherein the matrix elements in the resistivity variation matrix represent the resistivity ratio of the resistivity model at each grid unit at adjacent moments; According to the matrix value of the resistivity variation matrix, a first weight coefficient corresponding to each grid unit is determined.

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

[0052] Specifically, the resistivity variation matrix is ​​expressed as:

[0053] 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.

[0054] Furthermore, the step of determining the 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, the first weight coefficient corresponding to each grid unit is set 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.

[0055] For example, the first setting value may be set to 0.8, the second setting value may be set to 0.85, the third setting value may be set to 1.15, the fourth setting value may be set to 1.2, the first target value may be set to 0.01, and the second target value may be set to 1. For details, see Table 2:

[0056] Table 2. Correspondence between matrix value and first weight coefficient In an optional embodiment, 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.

[0057] In actual scenarios, since the grid setting range of the study area is larger than that of the actual observation area, the model change rate of the boundary grid cells can be ignored, so the second weight coefficients of the boundary grid cells are all 0.

[0058] Here, each non-boundary grid cell must obtain a variable result value according to the resistivity change relationship, and the variable result value of each non-boundary grid cell is different.

[0059] In an optional embodiment, the resistivity change relationship is expressed as:

[0060] 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.

[0061] In an optional embodiment, the step of 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 the 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, the second weight coefficient corresponding to each non-boundary grid unit is set to the fourth target value; If the variable result value is not greater than the second variable threshold, the second weight coefficient corresponding to each non-boundary grid unit is set 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.

[0062] Here, according to To determine the resistivity change of the resistivity model at each grid unit at the current moment, if The larger the value of is, the more likely it is that the resistivity model has changed at this location at the current moment. In this case, the grid cell is considered to be a gradual change area, and the second weight coefficient can be set to 1. The value of is medium, then the second weight coefficient can be set to 0.1; if If the value of is too small, it means that the resistivity model has no obvious change at this position at the current moment, then the grid unit is judged as a mutation area (whether there is a change at the next moment, it is not necessary to consider the change at the grid unit), and the second weight coefficient can be set to 0.01.

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

[0064] Table 3. Correspondence between variable result value and second weight coefficient In an optional embodiment, 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 characterizes the proportion of the model item in the objective function.

[0065] Here, the data term, model term and the second weight factor cooperate with each other in the time-lapse inversion and act together on the objective function to guide the optimization and adjustment of the resistivity model, and finally obtain a more accurate and reasonable underground model.

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

[0067] 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.

[0068] Here, the inversion is achieved by continuously minimizing the objective function, It will continue to change, and eventually the resistivity model that best fits the underground conditions will be obtained.

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

[0070] In this step, after multiple iterations, the objective function is minimized. When certain convergence conditions are met (such as the objective function value no longer decreases significantly, the parameter change is less than a certain threshold, etc.), the final resistivity model is obtained. The resistivity model describes the distribution of the resistivity of the underground medium. These models can intuitively display the electrical properties of different underground locations, providing important references for geological interpretation, resource exploration (such as finding mineral resources, groundwater resources, etc.), engineering construction (such as evaluating foundation stability, etc.), etc.

[0071] For example, in the experiment, at time T1, a small copper plate is selected as the abnormal body, that is, a small copper plate is inserted on the left side of the water tank; at time T2, a large copper plate is selected as the abnormal body, that is, a large copper plate is inserted on the left side of the water tank; at time T3, a combination of large and small copper plates is selected as the abnormal body, that is, a large copper plate and a small copper plate are inserted on the left and right sides of the water tank respectively. Under the above experimental conditions, we can get Figure 5 The inversion results are shown. Specifically, Figure 5The resistivity inversion result comparison diagram of the traditional method and the method of the present invention provided in the embodiment of the present application shows the resistivity inversion result profile at a depth of Z=0.15m, wherein white and black marking circles are used to mark the abnormal positions in the inversion results. It should be noted that since the black background cannot show the black marking circle, the white marking circle is used instead. It can be seen that the traditional method result (left) has anomalies at the left and right sides at three moments, and the position where the anomaly occurs at the third moment is significantly different from the actual position of the copper plate; the position where the anomaly occurs in the method of the present invention (right) is close to the actual position of the copper plate. In other words, when the traditional method is used for time-shift inversion, anomalies will appear on the right side of the first two moments, but there is actually no copper plate on the right side. At the same time, the position where the anomaly occurs at the third moment is significantly different from the actual position of the copper plate. The above-mentioned abnormalities are all caused by improper setting of the weight factor during the time-shift inversion process. It can be seen that the method of the present invention can avoid this problem and has obvious advantages.

[0072] The time-shift inversion method based on time-shift term constraints provided in the embodiment of the present application selects appropriate weight factors for the time-shift terms, specifically, selects appropriate first weight coefficients and second weight coefficients according to different changes in the underground area, so as to better reflect the actual changes of the underground model, reduce the error of the inversion results, and thus improve the accuracy of the inversion results.

[0073] Based on the same inventive concept, the embodiment of the present application also provides a time-shift inversion device based on time-shift term constraints corresponding to the time-shift inversion method based on time-shift term constraints. Since the principle of solving the problem by the device in the embodiment of the present application is similar to the above-mentioned time-shift inversion method based on time-shift term constraints in the embodiment of the present application, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be repeated.

[0074] See also Figure 6 , Figure 6 The schematic diagram of the structure of a time-shift inversion device based on time-shift term constraints provided in an embodiment of the present application is shown in FIG. Figure 6 As shown in , the time-shift inversion device 700 comprises: The data acquisition module 701 is used to acquire the apparent resistivity data of the underground medium; A grid division module 702 is used to divide the study area corresponding to the underground medium into a plurality of grid cells, wherein 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; A function construction module 703 is used to construct an objective function of time-shift inversion, wherein the objective function includes a time-shift term and a first weight factor corresponding to the time-shift term, wherein 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, wherein the first weight coefficient is used to reflect the resistivity change of the resistivity model at each grid unit at adjacent moments, and the second weight coefficient is 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 model updating module 704 is used to utilize the apparent resistivity data to continuously and iteratively update the model parameters of the resistivity model so as to continuously reduce the objective function and obtain an optimal resistivity model.

[0075] See also Figure 7 , Figure 7 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present application. Figure 7 As shown in FIG. 8 , the electronic device 800 includes a processor 801 , a memory 802 and a bus 803 .

[0076] The memory 802 stores machine-readable instructions executable by the processor 801. When the electronic device 800 is running, the processor 801 communicates with the memory 802 via the bus 803. When the machine-readable instructions are executed by the processor 801, the above-mentioned Figure 1 The specific implementation of the steps of the time-shift inversion method based on time-shift term constraints in the method embodiment shown can be found in the method embodiment, which will not be described in detail here.

[0077] The present application also provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the computer program can execute the above-mentioned Figure 1 The specific implementation of the steps of the time-shift inversion method based on time-shift term constraints in the method embodiment shown can be found in the method embodiment, which will not be described in detail here.

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

[0079] In the 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 schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation. For 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 mutual coupling or direct coupling or communication connection shown or discussed can be through some communication interfaces, and the indirect coupling or communication connection of devices or units can be electrical, mechanical or other forms.

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

[0081] In addition, each functional unit in each embodiment of the present application may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit.

[0082] If the functions are implemented in the form of software functional units and sold or used as independent products, they can be stored in a non-volatile computer-readable storage medium that is executable by a processor. Based on this understanding, the technical solution of the present application, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium, including several instructions for a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present application. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM), random access memory (RAM), disk or optical disk, and other media that can store program codes.

[0083] Finally, it should be noted that the above-described embodiments are only specific implementation methods of the present application, which are used to illustrate the technical solutions of the present application, rather than to limit them. The protection scope of the present application is not limited thereto. Although the present application is described in detail with reference to the above-mentioned embodiments, ordinary technicians in the field should understand that any technician familiar with the technical field can still modify the technical solutions recorded in the above-mentioned embodiments within the technical scope disclosed in the present application, or can easily think of changes, or make equivalent replacements for some of the technical features therein; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application, and should be included in the protection scope of the present application. Therefore, the protection scope of the present application shall be based on 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 cell at the current moment, represents the resistivity value of the resistivity model at the front grid unit at the current moment, Indicates 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, Indicates 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