Magnetotelluric probability inversion method based on image compression

Abstract

The invention relates to the technical field of electromagnetic exploration, in particular to a magnetotelluric probability inversion method based on image compression. According to the method, order reduction of the parameter model is achieved through image compression, so that the number of parameters is reduced, the number of the parameters is changed in a self-adaptive mode through the tower type structure, and the influence of human constraints on parameter updating in the inversion process is effectively avoided. Meanwhile, in the iteration process, the forward modeling speed is accelerated through process parallel, the inversion speed is accelerated through thread parallel, and therefore the sampling process of inversion is simpler, more convenient and more efficient. Finally, the magnetotelluric high-dimensional probability inversion is successfully realized, the inversion result can be quantitatively evaluated, and efficient and large-scale electromagnetic exploration can be realized.

Description

Technical Field

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

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

[0003] In electromagnetic inversion, due to the serious measurement interference, the inversion is easily affected by noise, which makes the inversion results unstable. 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 sense, and cannot reflect the non-uniqueness of 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 be able to directly extract models from different polarization data and perform quantitative evaluation of the inversion solution of the extracted model, a magnetotelluric probability inversion method based on image compression is urgently needed. 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 achieve the reduction of the order of the parameter model, thereby reducing the number of parameters, and uses a tower structure to adaptively change the number of parameters, effectively avoiding the influence of artificial constraints on parameter updates during the inversion process. At the same time, in the iterative process, the process parallelism is used to speed up the forward modeling speed, and the thread parallelism is used to speed up the inversion speed, so that the inversion sampling process is simpler and more efficient, so that the distribution of the posterior probability model related to the observation data can be efficiently obtained, and the quantitative evaluation of the magnetotelluric high-dimensional inversion can be realized. The specific technical scheme is as follows: A magnetotelluric probability inversion method based on image compression comprises the following steps: 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, and the tower structure of the initial sampling is set to start from only one node; 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 inverse transformation on the new candidate model under the tower structure obtained in S3, transform the inversion model in the wavelet domain into a suggested model of the conductivity structure, perform forward calculation on the suggested 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.

[0006] Optionally, 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 .

[0007] Optionally, in S2, the number of non-uniform tower permutations 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.

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

[0009] Optionally, in S3, the proposed model is transformed from the wavelet domain into a conductivity model in the physical space, the conductivity model including ,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.

[0010] Optionally, in S5, a 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. ; Optionally, in S5, a 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. .

[0011] The calculation of the proposed model acceptance rate is as follows: Select the proposed model acceptance rate of the existing node The expression is: ; Proposed 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.

[0012] Optionally, 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; If 0< <1, 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.

[0013] The application of the technical solution of the present invention has the following beneficial effects: (1) The present invention provides a method for magnetotelluric probability inversion based on image compression. The method of the present invention is based on the probability inversion of the Bayesian framework. It generates a set of experimental models through a probability sampling method, and performs statistical inference from the model set. It can reflect the non-uniqueness and uncertainty of the inversion results and make a quantitative evaluation of the inversion results. 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 realize multi-frequency parallel computing under different polarization modes, accelerates the forward calculation speed, and greatly reduces the time consumption of forward sampling of millions of levels.

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

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

[0016] (4) The method of the present invention avoids sampling from falling into local extreme values ​​through the parallel tempering method, and sets up thread parallel technology when multiple chains are tempered in parallel to further accelerate the inversion speed. Compared with the existing technology, the method of the present invention can effectively reduce the number of model inversions and accelerate the forward and inversion speed, thereby realizing the inversion using two-dimensional magnetotelluric different response data and quantitatively evaluating the inversion results.

[0017] In addition to the above-described objects, features and advantages, the present invention has other objects, features and advantages. The present invention will be further described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] In order to more clearly illustrate the embodiments of the present invention or the technical solutions of the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0019] Figure 1 is a flowchart of the steps of the magnetotelluric probability inversion method in a preferred embodiment of the present invention; Figure 2 is a schematic diagram of a non-uniform tower structure in a preferred embodiment of the present invention; Figure 3 is a schematic diagram of the exchange of multiple Markov chains in a preferred embodiment of the present invention; Figure 4 is a schematic diagram of an underground conductivity model in a preferred embodiment of the present invention, wherein triangles represent the locations of measuring stations; Figure 5 is a schematic diagram of a mean value model in a preferred embodiment of the present invention, wherein triangles represent the locations of measuring stations; FIG6 (a) and FIG6 (b) are one-dimensional edge probability density distribution diagrams obtained by sampling results at different measuring stations in a preferred embodiment of the present invention; FIG. 7 (a) is a schematic diagram showing the variation of the fitting error of the Markov chain at different temperatures with the number of iterations in a preferred embodiment of the present invention; FIG7( b ) is a schematic diagram showing the statistical distribution of the fitting errors of Markov chains at different temperatures in a preferred embodiment of the present invention; FIG8 (a) shows the variation of the number of nodes of the Markov chain at different temperatures with the number of iterations in the preferred embodiment of the present invention; FIG8( b ) is a probability distribution of the number of nodes of Markov chains at different temperatures in a preferred embodiment of the present invention; Fig. 9 is an actual inversion parameter diagram at a certain sampling moment during the sampling process in the preferred embodiment of the present invention; Fig.10 It is a comparison chart of the running time of the parallel structure running 10,000 times and the serial structure running 10,000 times in the preferred embodiment of the present invention. DETAILED DESCRIPTION

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

[0021] In order to enable those skilled in the art to better understand the scheme of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0022] like Figure 1 As shown, this embodiment provides a magnetotelluric probability inversion method based on image compression, comprising the following steps: 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 , maximum number of nodes Specifically, the tower structure in this embodiment is as follows: Figure 2 As shown, the maximum depth of the tower structure illustrated in this embodiment is 3.

[0023] 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. In this embodiment, the temperature gradient is set to 1.5, and the temperature of different Markov chains can be expressed as , n represents the total number of Markov chains. And multiple Markov chains are sampled simultaneously using thread parallel technology.

[0024] Furthermore, the number of non-uniform three-node tower arrangement combinations in this embodiment 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.

[0025] 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 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 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 tower structure in the wavelet domain, calculate the acceptance rate of the proposed model, and choose to accept or reject the new candidate model according to the acceptance rate; S6: Use the parallel tempering method to fully explore the space based on Markov chains at different temperatures; for each n samplings, randomly select two Markov chains at different temperatures, and determine whether to exchange the proposed models obtained by sampling the two 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; in this embodiment, 10≤Z≤20, Q is not less than 1000000, and Z and Q are natural numbers.

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

[0027] Specifically, 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.

[0028] Further, in S3, the proposed model is transformed from the wavelet domain into a conductivity model of the physical space, and the conductivity model expression is: ,in, represents the conductivity model after inverse transformation, represents the inverse wavelet transform operator; the proposed model is used for forward calculation. In this embodiment, in order to utilize all electromagnetic data forms, the proposed model The forward calculation includes the TE polarization mode and TM polarization mode , , and the method of this embodiment uses process parallel technology to reduce the time of obtaining magnetotelluric response data information obtained by different propagation characteristics of natural fields, that is, simultaneously obtaining magnetotelluric response data of different polarization modes at different frequencies, thereby reducing the forward calculation time of this step, wherein: 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.

[0029] It should be noted that in S4, the present embodiment uses MPI parallel technology to realize synchronous calculation of different processes. Specifically, the present embodiment uses MPI parallel technology to allocate different frequencies under different models to different processes for synchronous calculation, thereby accelerating the forward calculation time. The present embodiment can allocate the number of computer processes according to the number of frequency points that need to be calculated, such as the need to calculate the forward electromagnetic response of 4 frequencies under the TE polarization mode and the TM polarization mode. The present embodiment can provide 9 processes (process numbers are 0-8), and then process number 0 is regarded as the main process, which is used 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.

[0030] Optionally, in S5, a first acceptance rate of the proposed model is calculated ,like , then accept the new candidate model, that is, get the wavelet domain parameter model under the existing tower structure at the t+1th sampling Otherwise, the new candidate model is rejected, that is, the wavelet domain parameter model under the existing tower structure at the t+1th sampling is obtained. .

[0031] Optionally, in S5, a 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, which is preferably ; 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. .

[0032] The calculation of the proposed model acceptance rate is as follows: Select the proposed model acceptance rate of the existing node The expression is: ; Proposed 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.

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

[0034] Furthermore, in S6, in order to avoid falling into a local extreme value, the method of this embodiment exchanges the suggested models obtained by sampling the Markov chain at different temperatures, and the process is as follows: This embodiment uses OpenMP parallel technology to control Markov chains of different temperatures based on different threads and perform sampling at the same time, such as Figure 3 As shown, different temperatures Representing different threads, after n 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; If 0< <1, 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.

[0035] In order to verify the inversion effect of the method proposed in this embodiment, and thus prove that the required uncertainty information can be obtained, a model test is performed. Figure 4 The underground conductivity model shown in the figure, where the triangle represents the location of the measuring station, 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 , the conductivity of the abnormal body is The frequencies used are 8 frequency values: 0.1, 0.25, 1.0, 2.0, 5.0, 12.5, 20.0, 100 Hz.

[0036] The cross-dimensional sampler of this embodiment is implemented by programming in Fortran language, and is calculated using a high-performance platform, utilizing 9 cores. Figure 5This is the mean model of the underground medium distribution obtained by the cross-dimensional sampler of this embodiment, where the triangle represents the position of the measuring station. It can be seen that it is basically consistent with the actual model, and the size and position of the abnormal body can be well reflected. Figures 6 (a) and 6 (b) respectively represent the one-dimensional edge probability density distribution obtained by the sampling results of this embodiment at different measuring stations. The depth of the shadow color of the edge probability density distribution represents the conductivity probability of a certain area at the depth, and the darker the color, the higher the probability. From the results obtained by inversion, whether in the range of the abnormal body or in the background area, the results are consistent with the real 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. 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 static sampling is performed within the range, indicating that it has good mixing. Figure 7 (b) shows the statistical distribution of the fitting difference under different parallel tempering chains. As the temperature T increases, its distribution range becomes wider, indicating that its parallel tempering performance is better. 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, which shows that they have found the same inversion target, and only dozens of nodes can represent 1024 inversion parameters, indicating that it has good parameter compression. Figure 8 (b) shows the statistical distribution of the number of nodes under different parallel tempering chains. The distribution of different chains is basically the same, indicating that it has good convergence. Fig. 9 It can be seen that in the inversion process, the actual inversion node distribution is also obvious. It can also be clearly found that the inversion parameter model required for the original inversion is 32. 32=1024, only dozens of parameters are needed 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). Fig.10 It can be clearly seen that the parallel structure algorithm designed by the present invention is highly efficient and the calculation speed can be increased by 3.3 times.

[0037] The embodiment of the present invention provides a magnetotelluric probability inversion method based on image compression, which accelerates the forward modeling speed through the parallel tempering method, and uses the parallel tempering method to avoid sampling from falling into local extreme values. Compared with the prior art, the method of the present invention can effectively reduce the number of model inversions and speed up the forward and inversion speeds, thereby realizing the inversion using two-dimensional magnetotelluric different response data, and quantitatively evaluating the obtained inversion results.

[0038] It should be noted that the device embodiments described above are merely illustrative, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0039] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in 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 calculation of the proposed model acceptance rate is 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