A Magnetotelluric Probabilistic Inversion Method Based on Image Compression

Abstract

The present invention relates to the technical field of electromagnetic exploration, and specifically to a magnetotelluric probabilistic inversion method based on image compression. The method of the present invention uses image compression to reduce the order of the parameter model, thereby reducing the number of parameters, and adaptively changes the number of parameters using a tower structure, effectively avoiding the influence of artificial constraints on parameter updates during the inversion process. At the same time, during the iterative process, process parallelism is used to accelerate the forward modeling speed, and thread parallelism is used to accelerate the inversion speed, making the sampling process of the inversion more simple and efficient. Finally, the present invention successfully realizes magnetotelluric high-dimensional probabilistic inversion and can quantitatively evaluate the inversion results, which helps to achieve efficient and large-scale electromagnetic exploration.

Description

Technical Field

[0001] The present invention relates to the technical field of electromagnetic exploration, and particularly relates to a magnetotelluric probability inversion method based on image compression. Background Technique

[0002] Electromagnetic exploration plays an important role in mineral exploration, geological exploration, deep structure exploration, and underground metal exploration. Since the underground medium is non-uniform, two polarized electromagnetic response curves can be obtained at each measuring point of magnetotelluric sounding. One is the curve with the electric field polarized along the strike of the structure, which is called the TE-mode polarized electromagnetic response curve, and the other is the curve with the electric field polarized along the dip of the structure, which is called the TM-mode polarized electromagnetic response curve. The two polarized electromagnetic response curves can be used to describe different electromagnetic characteristics of the underground structure. Therefore, it has important practical application significance to invert the electromagnetic responses obtained under different polarization modes.

[0003] In electromagnetic inversion, due to relatively serious measurement interference, the inversion is easily affected by noise, resulting in unstable inversion results. There are multiple models fitting the measured data, and the inversion is non-unique. If a stable solution is obtained by forcibly adding artificial information, this solution is only the best-fitting model in a certain sense and cannot reflect the non-uniqueness characteristic existing in the inversion itself. Therefore, it is of great significance to combine different polarization modes and quantitatively evaluate the inversion solution only from the obtained response data.

[0004] In summary, in order to directly extract a model from different polarization data and quantitatively evaluate the inversion solution of the extracted model, there is an urgent need for a magnetotelluric probability inversion method based on image compression. Summary of the Invention

[0005] The purpose of the present invention is to provide a magnetotelluric probability inversion method based on image compression. The method of the present invention uses image compression to reduce the order of the parameter model, thereby reducing the number of parameters, and adaptively changes the number of parameters using a tower structure, effectively avoiding the influence of artificial constraints on parameter update during the inversion process. At the same time, during the iteration process, process parallelism is used to accelerate the forward modeling speed, and thread parallelism is used to accelerate the inversion speed, making the sampling process of the inversion simpler and more efficient, so as to efficiently obtain the posterior probability model distribution related to the observed data and realize the quantitative evaluation of magnetotelluric high-dimensional inversion. The specific technical solutions are as follows:

[0006] A magnetotelluric probability inversion method based on image compression, comprising the following steps:

[0007] S1: Receive the magnetotelluric response data information, construct a tower structure required for inversion in the wavelet domain, and calculate the maximum number of nodes of the tower structure based on the maximum depth of the tower structure;

[0008] S2: Calculate the tower permutation and combination numbers under nodes of different depths, set the temperature gradient, start constructing Markov chains at different temperatures from the first sampling, the Markov chains at different temperatures are parallel in threads, realize independent sampling from each other, and set the starting sampling tower structure to start from only one node;

[0009] S3: Divide the nodes of the tower structure into existing nodes, generable nodes, disposable nodes and empty nodes, select one type of node from the existing nodes, generable nodes and disposable nodes to change the tower structure, and form a new candidate model;

[0010] S4: Perform inverse wavelet transform on the new candidate model under the tower structure obtained in S3, transform the inversion model in the wavelet domain to obtain the proposed model of the conductivity structure, perform forward calculation on the proposed models corresponding to multiple Markov chains, and use the process parallel technology to obtain the magnetotelluric response data of different polarization modes at different frequencies, and accelerate the forward calculation process;

[0011] S5: Construct the Metropolis - Hasting acceptance criterion based on the tower structure in the wavelet domain, calculate the first acceptance rate of the proposed model, and select to accept or reject the new candidate model according to the first acceptance rate;

[0012] S6: After every n samplings, use the parallel tempering method, randomly select two Markov chains at different temperatures, and judge whether to exchange the proposed models obtained by sampling according to the second acceptance rate;

[0013] S7: Repeat steps S3 to S6, set the maximum number of samplings Q, save a new candidate model every Z samplings, and obtain multiple posterior probability models;

[0014] S8: Perform statistical calculations on all the posterior probability models obtained by sampling in S7, obtain all the information about the model, and complete the magnetotelluric inversion.

[0015] Optionally, in S1, the tower structure is a non - uniform tower structure, and the maximum depth of the non - uniform tower structure is , and the maximum number of nodes .

[0016] Optionally, in S2, the calculation formula of the non - uniform tower permutation and combination number is as follows:

[0017] ;

[0018] Among them, represents the depth of the tower structure, represents the number of tower nodes, represents at depth, The number of uniform tower - type structure arrangements under a node.

[0019] Optionally, in S3, the formation of the new candidate model includes:

[0020] Selecting the nodes that can generate for tower - type structure change, which represents the number of nodes after the change in the tower - type structure , the new candidate model is , where represents the number of nodes before the change, represents the wavelet coefficient value;

[0021] Selecting the nodes that can disappear for tower - type structure change, which represents the number of nodes in the new candidate tower - type structure , the new candidate model is ;

[0022] When selecting one of the existing nodes for tower - type structure change, it means only changing the value of the selected th node, while the number of nodes remains unchanged, , the new candidate model is , where represents the th new wavelet coefficient value of the node.

[0023] Optionally, in S3, the proposed model is transformed into a conductivity model in the physical space from the wavelet domain. The conductivity model includes , where represents the conductivity model after the inverse transformation, represents the wavelet inverse transformation operator; the proposed model is used for forward calculation, and the forward calculation of the proposed model includes the TE polarization mode and the TM polarization mode , , where:

[0024] The calculation process of the TE polarization mode is as follows:

[0025] ;

[0026] The calculation process of the TM polarization mode is as follows:

[0027] ;

[0028] Among them, is the electric field strength, is the magnetic field strength, is the medium conductivity, is the magnetic permeability of the medium, , is the frequency;

[0029] Calculate the likelihood ratio of the proposal model , and the expression is as follows:

[0030] ;

[0031] where represents the temperature factor.

[0032] Optionally, in S5, calculate the first acceptance rate of the proposal model , if , then accept the new candidate model, that is, obtain the existing tower structure at the (t + 1)-th sampling , otherwise reject the new candidate model, that is, obtain the existing tower structure at the (t + 1)-th sampling ;

[0033] Optionally, in S5, the process of calculating the first acceptance rate is as follows:

[0034] S5.1: Calculate the prior distribution of the current existing tower structure under the given tower structure and the probability distribution of the new candidate model ;

[0035] S5.2: Calculate the resistivity model obtained by inverse transformation of the existing tower structure model at a certain sampling time t likelihood ratio and the resistivity model obtained by inverse transformation of the new candidate tower structure model at a certain sampling time t likelihood ratio ;

[0036] S5.3: Calculate the proposal distribution probability from the existing tower structure model to the new candidate tower structure model and the proposal distribution from the new candidate tower structure model to the existing tower structure model ;

[0037] S5.4: Construct an acceptance criterion based on the prior distribution, likelihood ratio, and proposal distribution probability calculated in S5.1 to S5.3, and calculate the first acceptance rate of the proposal model. The expression of the first acceptance rate is as follows:

[0038] ;

[0039] Among them, represents the Jacobian matrix;

[0040] S5.5: According to the first acceptance rate calculated in S5.4 make a judgment. If 0 < < 1, then accept the proposed model, that is , otherwise reject the proposed model, that is .

[0041] The calculation of the acceptance rate of the proposed model is as follows:

[0042] Select the acceptance rate of the proposed model for the existing nodes The expression is:

[0043] ;

[0044] Select the acceptance rate of the proposed model for the producible nodes The expression is:

[0045] ;

[0046] Select the acceptance rate of the proposed model for the perishable nodes The expression is:

[0047] ;

[0048] Among them, represents the number of producible node sets of the current model, represents the number of perishable node sets of the current model, represents the number of producible node sets of the candidate model, represents the number of perishable node sets of the candidate model.

[0049] Optionally, in S6, exchange the proposed models obtained by Markov chain sampling at different temperatures. The process is as follows:

[0050] After multiple samplings, randomly select two Markov chains at different temperatures, and use the acceptance criterion to judge whether to exchange the proposed models under the two chains. The expression of the acceptance criterion is as follows:

[0051] ;

[0052] Among them, represents the second acceptance rate, represents the temperature of the proposed model, represents the temperature of the proposed model;

[0053] If 0 < < 1, then it is agreed that the two chains are exchanged, that is, the temperature under and the temperature under are exchanged, , otherwise the exchange is refused;

[0054] After the parallel tempering is completed, each Markov chain independently performs thread parallel sampling.

[0055] Applying the technical solution of the present invention has the following beneficial effects:

[0056] (1) The present invention provides a magnetotelluric probability inversion method based on image compression. The method of the present invention is based on the probability inversion of the Bayesian framework. By generating an experimental model set through a probability sampling method and performing statistical inference from the model set, the non-uniqueness and uncertainty of the inversion result can be reflected, and a quantitative evaluation can be made on the inversion result. Under the traditional Bayesian framework, there are problems such as long sampling time and a large number of magnetotelluric high-dimensional inversion parameters. The method of the present invention uses MPI parallel technology to achieve multi-frequency parallel calculation under different polarization modes, accelerating the forward calculation speed and greatly reducing the time-consuming of forward sampling of millions of samples.

[0057] (2) The method of the present invention uses image compression calculation to greatly reduce the number of inversion parameters, compressing thousands of inversion parameters into dozens, greatly reducing the degree of freedom of inversion.

[0058] (3) The method of the present invention uses a tower structure to adaptively change the number of model parameters according to the data complexity, making the inversion more flexible.

[0059] (4) The method of the present invention avoids sampling falling into local extrema through the parallel tempering method, and sets thread parallel technology during parallel tempering of multiple chains, further accelerating the inversion speed. Compared with the prior art, the method of the present invention can effectively reduce the number of model inversions, accelerate the forward and inverse inversion speeds, thereby realizing inversion using two-dimensional magnetotelluric different response data and making a quantitative evaluation of the obtained inversion result.

[0060] In addition to the purposes, features and advantages described above, the present invention has other purposes, features and advantages. The present invention will be further described in detail below with reference to the drawings. Brief Description of the Drawings

[0061] To more clearly illustrate the technical solutions of the embodiments of the present invention or the prior art, the following will briefly introduce the accompanying drawings required for the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other accompanying drawings can be obtained based on these drawings.

[0062] Figure 1 is the step flowchart of the magnetotelluric probability inversion method in the preferred embodiment of the present invention;

[0063] Figure 2 is the schematic diagram of the non-uniform tower structure in the preferred embodiment of the present invention;

[0064] Figure 3 is the schematic diagram of the exchange of multiple Markov chains in the preferred embodiment of the present invention;

[0065] Figure 4 is the schematic diagram of the underground conductivity model in the preferred embodiment of the present invention, where the triangles represent the positions of the measuring stations;

[0066] Figure 5 is the schematic diagram of the mean model in the preferred embodiment of the present invention, where the triangles represent the positions of the measuring stations;

[0067] Figures 6(a) and 6(b) are the one-dimensional marginal probability density distribution diagrams obtained from the sampling results at different measuring stations in the preferred embodiment of the present invention;

[0068] Figure 7(a) is the schematic diagram of the change of the fitting error of the Markov chains at different temperatures with the number of iterations in the preferred embodiment of the present invention;

[0069] Figure 7(b) is the schematic diagram of the statistical distribution of the fitting errors of the Markov chains at different temperatures in the preferred embodiment of the present invention;

[0070] Figure 8(a) is the change of the number of nodes of the Markov chains at different temperatures with the number of iterations in the preferred embodiment of the present invention;

[0071] Figure 8(b) is the probability distribution of the number of nodes of the Markov chains at different temperatures in the preferred embodiment of the present invention;

[0072] Figure 9 is the actual inversion parameter diagram at a certain sampling moment during the sampling process in the preferred embodiment of the present invention;

[0073] Figure 10 is the comparison diagram of the running time of the parallel structure running 10,000 times and the serial running 10,000 times in the preferred embodiment of the present invention. Detailed implementation manners

[0074] In high-dimensional magnetotelluric probabilistic inversion, due to the large number of parameters, it is difficult to sample the posterior probability density distribution of the model in high-dimensional space. To achieve high-dimensional magnetotelluric probabilistic inversion, the method of the present invention is based on the idea of image compression and adopts wavelet image compression technology to transform the inversion space from the physical space model into the wavelet domain model space, so that the inversion parameters have different scale characteristics and the number of inversion parameters is reduced. Combined with the tower structure, through the electromagnetic response data obtained based on different propagation characteristics, the self-feedback technology is used to intelligently determine the parameter changes, thereby realizing high-dimensional magnetotelluric inversion. In addition, to further improve the efficiency, the method of the present invention also adopts a process parallel technology to reduce the time for obtaining the magnetotelluric response data information of the natural field under different propagation characteristics. In the inversion process, the present invention introduces a parallel tempering technology to avoid sampling falling into local extreme points and accelerates the calculation process through thread parallelism.

[0075] In order to enable those skilled in the art to better understand the solution of the present invention, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without making creative efforts shall fall within the protection scope of the present invention.

[0076] As Figure 1 shown, this embodiment provides a magnetotelluric probabilistic inversion method based on image compression, including the following steps:

[0077] S1: Receive the magnetotelluric response data information, construct the tower structure required for inversion, and calculate the maximum number of nodes of the tower structure based on the maximum depth of the tower structure; specifically, the tower structure in this embodiment is a non-uniform tower structure, and the maximum depth of the non-uniform tower structure is , and the maximum number of nodes . Specifically, the tower structure in this embodiment is as Figure 2 shown, and the maximum depth of the tower structure exemplified in this embodiment is 3.

[0078] S2: Calculate the tower permutation and combination numbers at nodes with different depths, set the temperature gradient, and construct Markov chains at different temperatures starting from the first sampling; the temperature gradient in this embodiment is set to 1.5, and the temperatures of different Markov chains can be expressed as , where n represents the total number of Markov chains. And multiple Markov chains use thread parallel technology for simultaneous sampling.

[0079] Furthermore, the calculation formula for the non-uniform three-node tower permutation and combination number in this embodiment is as follows:

[0080] ;

[0081] Among them, represents the depth of the tower structure, represents the number of tower nodes, represents at depth, the number of uniform tower structure arrangements under nodes.

[0082] S3: Divide the nodes of the tower structure into existing nodes, generable nodes, perishable nodes, and empty nodes, select one type of node among the existing nodes, generable nodes, and perishable nodes for tower structure changes to form a new candidate model;

[0083] S4: Perform inverse wavelet transform on the new candidate model under the tower structure obtained in S3, convert the inversion model in the wavelet domain to obtain a proposed model of the conductivity structure, perform forward calculation on the proposed models corresponding to multiple Markov chains, and use process parallel technology to obtain magnetotelluric response data of different polarization modes at different frequencies to accelerate the forward calculation process;

[0084] S5: Construct a Metropolis - Hasting acceptance criterion based on the tower structure in the wavelet domain, calculate the acceptance rate of the proposed model, and select to accept or reject the new candidate model according to the acceptance rate;

[0085] S6: Use the parallel tempering method to fully explore the space based on Markov chains at different temperatures; for each n - th sampling, randomly select two Markov chains at different temperatures, and judge whether to exchange the proposed models obtained by sampling them according to the second acceptance rate;

[0086] S7: Repeat steps S3 to S6, set the maximum number of samplings Q, save a new candidate model every Z samplings, and obtain multiple posterior probability models; in this embodiment, 10 ≤ Z ≤ 20, Q is not less than 1000000, and Z and Q are natural numbers.

[0087] S8: Perform statistical calculations on all the posterior probability models obtained by sampling in S7 to obtain all information about the model and complete the magnetotelluric inversion.

[0088] Specifically, in S3, the formation of the new candidate model includes:

[0089] Select generable nodes for tower structure changes, then represents the number of changed nodes in the tower structure , the new candidate model is , among which, represents the number of nodes before the change, represents the wavelet coefficient value;

[0090] Selecting a perishable node for tower structure change indicates the number of new candidate tower structure nodes , and the new candidate model is ;

[0091] When selecting one of the existing nodes for tower structure change, it means only changing the value of the selected th node, while the number of nodes remains unchanged, , and the new candidate model is , where represents the new wavelet coefficient value of the th node.

[0092] Furthermore, in S3, the proposed model is transformed from the wavelet domain into a conductivity model in physical space, and the expression of the conductivity model is: , where represents the conductivity model after inverse transformation, represents the wavelet inverse transform operator; the proposed model is used for forward calculation. In this embodiment, in order to be able to utilize all data forms of electromagnetics, the forward calculation of the proposed model includes the TE polarization mode and the TM polarization mode , , and the method of this embodiment utilizes the process parallel technology to reduce the time for obtaining the magnetotelluric response data information obtained from different propagation characteristics of the natural field, that is, simultaneously obtaining the magnetotelluric response data of different polarization modes at different frequencies, thereby reducing the forward calculation time of this step, where:

[0093] The calculation process of the TE polarization mode is as follows:

[0094] ;

[0095] The calculation process of the TM polarization mode is as follows:

[0096] ;

[0097] where is the electric field strength, is the magnetic field strength, is the medium conductivity, is the magnetic permeability of the medium, , is the frequency;

[0098] Calculating the likelihood ratio of the proposed model , and the expression is as follows:

[0099] ;

[0100] Among them, represents the temperature factor.

[0101] It should be noted that in S4, the MPI parallel technology is adopted in this embodiment to realize the synchronous calculation of different processes. Specifically, in this embodiment, the MPI parallel technology is used to allocate different frequencies under different models to different processes for synchronous calculation, so as to accelerate the forward calculation time. This embodiment can allocate the number of processes of the computer according to the number of frequency points to be calculated. For example, it is necessary to calculate the forward electromagnetic responses at 4 frequencies in the TE polarization mode and the TM polarization mode. This embodiment can give 9 processes (process numbers are 0-8), and then the 0th process is regarded as the main process to summarize information, and the remaining process numbers are numbered according to odd and even numbers to calculate the TE polarization mode and the TM polarization mode respectively.

[0102] Optionally, in S5, calculate the first acceptance rate of the proposed model , if , then accept the new candidate model, that is, obtain the wavelet domain parameter model under the existing tower structure at the (t + 1)-th sampling , otherwise reject the new candidate model, that is, obtain the wavelet domain parameter model under the existing tower structure at the (t + 1)-th sampling .

[0103] Optionally, in S5, the process of calculating the first acceptance rate is as follows:

[0104] S5.1: Calculate the prior distribution of the current existing tower structure under the given tower structure and the probability distribution of the new candidate model

[0105] ; S5.2: Calculate the likelihood ratio of the resistivity model obtained by inverse transformation of the existing tower structure model at a certain sampling time t and the likelihood ratio of the resistivity model

[0106] obtained by inverse transformation of the new candidate tower structure model at a certain sampling time t ; S5.3: Calculate the proposal distribution probability from the existing tower structure model to the new candidate tower structure model ;

[0107] S5.4: Construct an acceptance criterion based on the prior distribution, likelihood ratio, and proposed distribution probability calculated in S5.1 to S5.3, and calculate the first acceptance rate of the proposed model. , the first acceptance rate has the following expression:

[0108] ;

[0109] where, represents the Jacobian matrix, and preferably in this embodiment ;

[0110] S5.5: Make a judgment based on the first acceptance rate calculated in S5.4. If 0 < < 1, then accept the proposed model, that is , otherwise reject the proposed model, that is .

[0111] The calculation of the proposed model acceptance rate is as follows:

[0112] The expression for the proposed model acceptance rate of the selected existing node is:

[0113] ;

[0114] The expression for the proposed model acceptance rate of the selected producible node is:

[0115] ;

[0116] The expression for the proposed model acceptance rate of the selected perishable node is:

[0117] ;

[0118] where, represents the number of producible node sets of the current model, represents the number of perishable node sets of the current model, represents the number of producible node sets of the candidate model, represents the number of perishable node sets of the candidate model.

[0119] It should be noted that since the magnetotelluric inversion problem is a non - linear problem and there are multiple local extreme points in its solution space, therefore, using the parallel tempering method, different temperature chains can be used to fully explore the space.

[0120] Further, in S6, to avoid falling into local extrema, the method of this embodiment exchanges the proposed models obtained by Markov chain sampling at different temperatures. The process is as follows:

[0121] This embodiment adopts the OpenMP parallel technology to control the Markov chains at different temperatures based on different threads and perform sampling simultaneously. As Figure 3 shown, different temperatures represent different threads. After n samplings, randomly select two Markov chains at different temperatures and use the acceptance criterion to judge whether to exchange the proposed models under the two chains. The expression of the acceptance criterion is as follows:

[0122] ;

[0123] where represents the second acceptance rate, represents the temperature the proposed model under it, represents the temperature the proposed model under it;

[0124] If 0 < < 1, it is agreed to exchange the two chains, that is, the under temperature and the under temperature are exchanged, , otherwise the exchange is rejected;

[0125] After performing parallel tempering, each Markov chain independently performs thread-parallel sampling.

[0126] To verify the inversion effect of the method proposed in this embodiment and thus prove that the required uncertainty information can be obtained, model testing is carried out. Based on the above method, this embodiment designs an underground conductivity model as Figure 4 shown, where the triangles represent the measurement station positions. A cross-dimensional sampler is constructed based on the underground conductivity model, and the posterior probability model is obtained by sampling through the cross-dimensional sampler. The details are as follows: The detection range of the underground conductivity model is , the background conductivity is , and the conductivity of the anomaly is . The frequencies used are 8 frequency values: 0.1, 0.25, 1.0, 2.0, 5.0, 12.5, 20.0, 100 Hz.

[0127] The cross-dimensional sampler of this embodiment is implemented by programming in Fortran language, and calculations are carried out using a high-performance platform, using 9 cores. Figure 5The mean model of the underground medium distribution sampled by the cross-dimensional sampler of this embodiment, where the triangles represent the station positions. It can be seen that it is basically consistent with the actual model, and the size and position of the anomaly can be well reflected. Figures 6(a) and 6(b) respectively show the one-dimensional marginal probability density distributions obtained at different stations for the sampling results of this embodiment. The depth of the shadow color of the marginal probability density distribution represents the probability of the conductivity of a certain area at that depth, and the darker the color, the higher the probability. From the inversion results, whether within the anomaly range or in the background area, the results are consistent with the true distribution results. And through Figures 7(a), 7(b), 8(a) and 8(b), it can be seen that the sampling of the cross-dimensional sampler in this embodiment has good mixing and constraint on the number of parameters, and greatly reduces the number of parameters. Among them, Figure 7(a) shows that as the number of inversion iterations increases, the data fitting error decreases rapidly, and multiple parallel tempering chains basically drop to the same range and perform static sampling within the range, indicating good mixing. Figure 7(b) shows the statistical distribution of the fitting errors under different parallel tempering chains. As the temperature T increases, the distribution range becomes wider, indicating better performance of the parallel tempering. Figure 8(a) shows the change of the number of nodes under different parallel tempering chains with the number of iterations. The nodes of different tempering chains are basically in the same range, indicating that they converge to the same inversion target, and only dozens of nodes can represent 1024 inversion parameters, indicating good parameter compressibility. Figure 8(b) shows the statistical distribution of the number of nodes under different parallel tempering chains. The distributions of different chains are basically the same, indicating good convergence. Through Figure 9 It can be seen that during the inversion process, the actual distribution of the inversion nodes, and it can also be clearly found that the original inversion parameter model required for inversion is 32 32 = 1024, and only dozens of parameters are required in the wavelet domain under image compression (the gray in the figure represents its value, the darker the color, the larger the value, and the white part represents no parameters). Figure 10 It can be clearly seen that the parallel structure algorithm designed by the present invention has high efficiency, and the calculation speed can be increased by 3.3 times.

[0128] The embodiment of the present invention provides a magnetotelluric probability inversion method based on image compression. By using the parallel tempering method, the forward modeling speed is accelerated, and the parallel tempering method is used to avoid the sampling falling into local extrema. Compared with the prior art, the method of the present invention can effectively reduce the number of model inversions, speed up the forward and inverse modeling speeds, so as to realize the inversion using different response data of two-dimensional magnetotellurics and quantitatively evaluate the obtained inversion results.

[0129] It should be noted that the device embodiments described above are merely illustrative. 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 to multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0130] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, the present invention can have various changes and modifications. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A magnetotelluric probability inversion method based on image compression, characterized in that: The following steps are involved: S1: Receive the magnetotelluric response data information, construct the tower structure required for inversion in the wavelet domain, and calculate the maximum number of nodes of the tower structure based on the maximum depth of the tower structure; S2: Calculate the number of tower combinations at different depth nodes, set the temperature gradient, and build Markov chains at different temperatures starting from the first sampling. The Markov chains at different temperatures are threaded in parallel to achieve independent sampling. S3: Divide the nodes of the tower structure into existing nodes, generateable nodes, extinct nodes and empty nodes, and select one of the existing nodes, generateable nodes and extinct nodes to change the tower structure to form a new candidate model; S4: Perform wavelet inversion on the new candidate model under the tower structure obtained in S3, convert the inversion model in the wavelet domain to obtain the recommended model of the conductivity structure, perform forward calculation on the recommended models corresponding to multiple Markov chains, and use process parallel technology to obtain magnetotelluric response data of different polarization modes at different frequencies to accelerate the forward calculation process; S5: construct the Metropolis-Hasting acceptance criterion based on the wavelet domain tower structure, calculate the first acceptance rate of the proposed model, and choose to accept or reject the new candidate model according to the first acceptance rate; S6: After every n samplings, use the parallel tempering method to select two Markov chains at different temperatures, and determine whether to exchange the sampled suggestion models according to the second acceptance rate; S7: Repeat steps S3 to S6, set the maximum number of sampling times Q, save a new candidate model every Z sampling times, and obtain multiple posterior probability models; S8: Perform statistical calculations on the posterior probability models obtained from all samples in S7 to obtain all information about the model and complete the magnetotelluric inversion.

2. The magnetotelluric probability inversion method according to claim 1, characterized in that: In S1, the tower structure is a non-uniform tower structure, and the maximum depth of the non-uniform tower structure is , maximum number of nodes .

3. The magnetotelluric probability inversion method according to claim 2, characterized in that: In S2, the number of non-uniform tower permutations is The calculation formula is as follows: ； in, Expressed as the tower structure depth, Expressed as the number of tower nodes, Indicated in depth, The number of uniform tower structures under each node.

4. The magnetotelluric probability inversion method according to claim 3, characterized in that: In S3, the formation of a new candidate model includes: Select the node that can be generated to change the tower structure, which means the number of nodes in the tower structure after the change , the new candidate model is ,in, Indicates the number of nodes before the change, Expressed as wavelet coefficient value; Select the nodes that can be eliminated to change the tower structure, which means the number of new candidate tower structure nodes , the new candidate model is ; When one of the existing nodes is selected for tower structure change, it means that only the selected node is changed. The value of the nodes, while the number of nodes remains unchanged, , the new candidate model is ,in, Indicates The new wavelet coefficient value of each node.

5. The magnetotelluric probability inversion method according to claim 4, characterized in that: In S4, the proposed model is a conductivity model converted from the wavelet domain to the physical space, which includes ,in, represents the conductivity model after inverse transformation, represents the inverse wavelet transform operator, the proposed model is used for forward calculation, the proposed model The forward calculation includes the TE polarization mode and TM polarization mode , ,in: TE polarization mode The calculation process is as follows: ； TM polarization mode The calculation process is as follows: ； in, is the electric field strength, is the magnetic field strength, is the medium conductivity, is the magnetic permeability of the medium, , is the frequency; Calculate the likelihood ratio of the proposed model , the expression is as follows: ； in, Represents the temperature factor.

6. The magnetotelluric probability inversion method according to claim 5, characterized in that: In S5, the first acceptance rate of the proposed model is calculated ,like , then accept the new candidate model, that is, get the existing tower structure at the t+1th sampling Otherwise, the new candidate model is rejected, that is, the existing tower structure at the t+1th sampling time is obtained. .

7. The magnetotelluric probability inversion method according to claim 6, characterized in that: In S5, the first acceptance rate is calculated The process is as follows: S5.1: Compute the prior distribution of the current existing tower structure given a tower structure And the probability distribution of the new candidate model ; S5.2: Calculate the existing tower structure model at a certain sampling time t Resistivity model obtained by inverse transformation The likelihood ratio And the new candidate tower structure model at a certain sampling time t Resistivity model obtained after inverse transformation The likelihood ratio ; S5.3: Calculate the existing tower structure model To the new candidate tower structure model The proposed distribution probability And re-select tower structure models To the existing tower structure model Suggested distribution of ; S5.4: Construct acceptance criteria based on the prior distribution, likelihood ratio, and proposed distribution probability calculated in S5.1 to S5.3, and calculate the first acceptance rate of the proposed model , the first acceptance rate The expression is as follows: ； in, represents the Jacobian matrix; S5.5: The first acceptance rate calculated according to S5.4 Make a judgment, if 0< <1, then accept the proposed model, that is, , otherwise the proposed model is rejected, i.e. ; The proposed model acceptance rate is calculated as follows: Select the proposed model acceptance rate of the existing node The expression is: ； Proposal model acceptance rate for selecting production-ready nodes The expression is: ； Acceptance rate of the proposed model for selecting extinct nodes The expression is: ； in, Indicates the number of nodes that can be generated in the current model. Represents the number of extinct nodes in the current model. Represented as the number of generative node sets for candidate models, The number of nodes that can be eliminated as candidate models.

8. The magnetotelluric probability inversion method according to claim 7, characterized in that: In S6, the proposed models obtained by sampling the Markov chains at different temperatures are exchanged, and the process is as follows: After multiple samplings, two Markov chains at different temperatures are randomly selected, and the acceptance criterion is used to determine whether the proposed models under the two chains are exchanged. The expression of the acceptance criterion is as follows: ； in, represents the second acceptance rate, Indicates temperature The proposed model below, Indicates temperature The proposed model below; like , then the two chains are allowed to exchange, that is, the temperature Next With temperature Next exchange, , otherwise the exchange is refused; After executing parallel tempering, each Markov chain performs thread-parallel sampling independently.

Images