CN117574790A - Design method of cross-dimension Bayesian sampler based on physical space tree structure - Google Patents

Design method of cross-dimension Bayesian sampler based on physical space tree structure Download PDF

Info

Publication number
CN117574790A
CN117574790A CN202410078165.4A CN202410078165A CN117574790A CN 117574790 A CN117574790 A CN 117574790A CN 202410078165 A CN202410078165 A CN 202410078165A CN 117574790 A CN117574790 A CN 117574790A
Authority
CN
China
Prior art keywords
model
tree structure
node
nodes
distribution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202410078165.4A
Other languages
Chinese (zh)
Other versions
CN117574790B (en
Inventor
田圣琦
郭荣文
柳建新
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202410078165.4A priority Critical patent/CN117574790B/en
Publication of CN117574790A publication Critical patent/CN117574790A/en
Application granted granted Critical
Publication of CN117574790B publication Critical patent/CN117574790B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Physics (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Pure & Applied Mathematics (AREA)
  • Geometry (AREA)
  • Medical Informatics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Linguistics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a design method of a cross-dimensional Bayesian sampler based on a physical space tree structure, which relates to the technical field of geophysical inversion and comprises the following steps: measuring the response of the underground electromagnetic field of a certain area; constructing a tree structure initial model; optional tree structurem t Tree structure change is carried out on the class of nodes of the model to form a new candidate modelm'The method comprises the steps of carrying out a first treatment on the surface of the New candidate modelm'Inversely transforming to obtain a proposal model, and calculating forward modeling response; calculating a likelihood function; calculating acceptance rate of suggestion modelαJudgingαWhether or not it is smaller than a random numberIf it is, then there arem t+1 = m'If otherwise there ism t+1 =m t The method comprises the steps of carrying out a first treatment on the surface of the Taking outt=t+1Repeated sampling, per runZSubsampling holds the one-time sampling model untilt=QEnding the time; each will beZAnd carrying out statistical calculation on the sampling model stored by sub-sampling to obtain one-dimensional edge probability density distribution. The invention solves the problem that the prior distribution of the nodes of the wavelet domain tree structure is difficult to select, and improves the sampling efficiency.

Description

Design method of cross-dimension Bayesian sampler based on physical space tree structure
Technical Field
The invention relates to the technical field of geophysical inversion, in particular to a design method of a cross-dimensional Bayesian sampler based on a physical space tree structure.
Background
Geophysical magnetotelluric exploration is a geophysical method for acquiring underground physical response by using a natural electromagnetic field source, and is widely applied to research in the fields of oil gas, mineral products, underground deep-layer ground electrical structures and the like. The inversion results also show non-uniqueness due to interference in the actual observations. Therefore, along with the deep exploration and the rapid improvement of the development level of computer software and hardware, the reliability evaluation of the magnetotelluric inversion result gradually becomes the focus and hot spot of the study of domestic and foreign students. Uncertainty analysis is carried out on the inversion result, the distribution situation of possible underground electrical property is given, quantitative evaluation of the inversion result by people is greatly improved, and guidance is given in actual engineering.
At present, probability inversion based on a Bayesian framework is performed in order to embody uncertainty of an inversion result and quantitatively evaluate the inversion result. Model extraction is carried out by using a Markov Monte Carlo (Markov chain MonteCarlo, mcMC) sampling method, and statistical inference of parameters is carried out.
In bayesian inversion, the dimensions of model parameters cannot be determined in advance because the knowledge of the subsurface structure for actual exploration is limited. The literature (Green, p.j., 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, biomericka, 82 (4), 711-732) considers model parameters as unknowns and incorporates them into a bayesian framework, a method known as trans-dimensional bayesian inversion. The complexity of the model is driven by data through the cross-dimensional Bayesian inversion, and the optimization parameters do not need to be selected manually, so that the influence of human factors on the inversion result is avoided, and the method becomes a mainstream method. The literature (Mallingerno, A., 2002, parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem, geophys. J. Int., 151 (3), 675-688.) first discloses a natural extension of Bayesian parameter estimation in inverting the DC resistance, i.e., parameterizing the earth model into a layered medium, inverting the number of layers as parameters, and obtaining samples of the layered medium from a posterior distribution. By means of this automatic parameterization, it is possible to provide a higher model resolution in areas with better data coverage and low noise. However, this parameterization method cannot address the high-dimensional magnetotelluric inversion problem.
Currently, the Voronoi cell method is a common method among methods for performing parameter cross-dimensional sampling in a high-dimensional parameter space. Literature (Bodin, t., sambridge, m., rawlinson, n. & aroucau, p., 2012, transdimensional tomography with unknown data noise, geophys, j. Int., 189, 1536-1556.) discloses the discretization of a subsurface conductivity profile into a series of Voronoi polygons by constructing a set of continuous polygons by connecting perpendicular bisectors of adjacent nuclei. The cross-dimensional sampling of the number of model parameters is achieved by increasing or decreasing the number of nuclei. The method has the advantage that the parameter position and physical parameters of the geophysical model space can be represented by specifying the position of atomic nuclei and the attribute value inside each cell, so that the model is flexibly divided. However, the Voronoi cell method requires three unknowns to be contained in one cell. As the data volume of the geophysical model increases, the complexity and time cost of the geophysical model increases. Furthermore, the mesh subdivision of the Voronoi cell subdivision method is a non-uniform arbitrary subdivision, which presents challenges for the finite difference-based geoelectromagnetic forward computing method. Literature (Hawkins, r. & Sambridge, m. & 2015, geophysical imaging using transdimensional trees, geophys, j. Int, 203, 972-1000.) proposes a method for sampling wavelet tree structures. The method converts inversion of the geophysical model in physical space into research and operation in a tree structure based on wavelet basis functions, and can represent a complex geophysical model with a small number of parameters. However, since the wavelet coefficients are stored in the tree-like structure nodes and do not represent any physical significance, the method is difficult to choose prior distribution and cannot correspond to the conductivity range in the physical space. Meanwhile, the selection of different wavelet basis functions can directly influence the final inversion result, so that the application condition of the method is limited to a certain extent. Therefore, it is necessary to provide a new tree structure to solve the above technical problems.
Disclosure of Invention
The invention mainly aims to provide a design method of a cross-dimensional Bayesian sampler based on a physical space tree structure, which aims to solve the problems that prior distribution of nodes of a wavelet domain tree structure is difficult to select and sampling efficiency is low.
In order to achieve the above object, the present invention provides a design method of a cross-dimensional bayesian sampler based on a physical space tree structure, which includes:
s1: measuring by using an electromagnetic measuring instrument to obtain the response of an underground electromagnetic field of a certain area, and determining the measured underground depth;
s2: constructing a tree structure initial model, and constructing a Markov chain, wherein the method specifically comprises the following steps of:
calculating the maximum depth of the tree structure, wherein the formula is as follows:wherein->Representing the number of subdivisions of the subsurface depth at the maximum depth of the tree-like structure,/->Representing the maximum depth of the tree structure and constructing the maximum depth +.>A lower binary tree structure;
number of times of construction samplingTime tree structure initial model->Wherein->For the number of nodes in the initial model, +.>Depth level for tree structure, ++>， />Expressed as depth->Lower part (C)The number of arrangement of tree structures of all different shapes composed of the number of nodes, < >>Indicate->Conductivity parameters formed by the individual nodes;
from the slaveWhen starting to build Markov chain for +.>Existing tree structure at subsampling +.>In other words, the nodes of the tree structure can be classified into existing nodes, generatable nodes, and deceased nodes;
s3: optional tree structureTree structure change is carried out on the class of nodes of the model (1) to form a new candidate model +.>The method comprises the steps of carrying out a first treatment on the surface of the If a generatable node is selected, a new candidate tree node is represented + ->The new candidate model isThe method comprises the steps of carrying out a first treatment on the surface of the If an erasable node is selected, a new candidate tree node is indicated +.>The new candidate model is +.>The method comprises the steps of carrying out a first treatment on the surface of the If one of the existing nodes is selected, it means that only the conductivity parameter of the selected node is changed +.>The new candidate model is +.>；
S4: new candidate models under the tree structure obtained in the step S3Performing inverse transformation to obtain a proposed model of the conductivity structure>And by forward algorithm ∈ ->Will suggest a model->Mapping to a data space to obtain forward response;
s5: calculation based on the electromagnetic field response in S1 and the forward response obtained in S4Conference modelLikelihood functions of (2);
s6: obtaining acceptance rate of suggestion model based on likelihood function calculation in S5Taking random number between 0 and 1, and judging the acceptance rate +.>Whether the random number is smaller than the random number, if yes, receiving a new candidate model to obtain the +.>The existing tree structure at sub-sampling is: />Otherwise, rejecting the new candidate model to obtain +.>The existing tree structure at sub-sampling is:；
s7: order theRepeating S3 to S6 and performing +/every time>Sampling model obtained by sampling once is stored until sampling times +.>Ending the time;
s8: and (3) carrying out statistical calculation on all the sampling models obtained in the step (S7) to obtain the one-dimensional edge probability density distribution of the underground conductivity of the area.
Optionally, the inverse transformation in S4 specifically includes:
construction of a Scale functionTo control the currently existing nodes, the specific formula is as follows:
；
and for tree structure depth isLower->The individual nodes are defined as:
；
wherein,expressed as +.>Unit function on->Any number expressed as 0-1 has no specific meaning, < >>Expressed as the level at which the node is located, +.>The position of the level of the node is shown;
constructing a decoupling functionThe relationship between the upper level node and the next level node is decoupled, namely the value of the parent node of the range where the child node is located is covered by the value of the child node, and the specific formula is as follows:
；
similarly, for tree structure depth isLower->The individual nodes define:
；
wherein,expressed as +.>Unit function on the table.
Optionally, the S5 includes:
s5.1: suggested model obtained according to inverse transformation of tree-structured candidate modelLeast square method is carried out on the forward calculation response of (2) and the electromagnetic field response obtained by measurement, and fitting difference +.>The specific formula is as follows:
；
wherein,representing the mapping of the model vector to the positive operator in the data space, < >>Is the data covariance matrix,>response data for the actual measured electromagnetic field;
s5.2: assuming that likelihood functions obey a standard normal distribution, likelihood functionsThe specific formula of (2) is as follows:
；
wherein,for the number of observations.
Optionally, the S6 includes:
s6.1: computing a priori distribution of a current existing tree structure given a binary tree structurePrior distribution of new candidate models +.>；
S6.2: calculate the firstExisting tree-like structure model at subsampling +.>Inverse transformation of the resulting conductivity model->Likelihood function of +.>First->New candidate tree-like structure model at subsampling +.>Suggested model obtained by inverse transformationLikelihood function of +.>；
S6.3: calculation from existing tree structure modelTo a new candidate tree-structure model->Is a suggested distribution probability of (2)From the new candidate tree structure model +.>To the existing tree structure model->Is->；
S6.4: calculating the acceptance rate of the proposal model according to the prior distribution, likelihood function and proposal distribution calculated in S6.1 to S6.3The specific formula is as follows:
；
wherein,is Jacobian matrix->；
S6.5: taking a random number between 0 and 1, and judging the acceptance rateWhether or not it is smaller than the random number, if soThen accept the advice model, i.e->Otherwise reject the advice model, i.e. +.>。
Optionally, the step S6.1 includes:
s6.1.1: calculating uniform distribution over a given range of node numbersThe specific formula is as follows:
；
wherein,is the maximum node number of the tree structure, +.>；/>Is the minimum node number of the tree structure, +.>Generally set to 1, indicating at least one node is owned;
s6.1.2: calculating the prior number of the arrangement under the given node numberNamely, the distribution of the permutation and combination number of the tree structure under the condition of a given node is as follows:
；
wherein,is->The number of nodes forms the number of tree-like structures of all different shapes,，/>is a standard binomial coefficient;
s6.1.3: computing a priori distribution of parameters in each given nodeThe specific formula is as follows:
；
s6.1.4: the prior distribution under the binary tree structure is calculated respectively by adopting the same formulaAnd->Wherein a priori distribution->The specific formula of (2) is as follows:
。
optionally, the S6.3 includes:
s6.3.1: forming the current nodes on the current tree structure into a perturbable node groupThe generatable nodes form a generatable node group +.>The eroding nodes form eroding node group +.>；
S6.3.2: if the node selected in the S3 is the existing node, calculating disturbance suggestion distribution, if the node selected in the S3 is a generatable node, calculating generation suggestion distribution, and if the node selected in the S3 is an deceaseable node, calculating deceased suggestion distribution.
Optionally, the calculating of the disturbance recommendation distribution in S6.3.2 includes:
calculating a givenSelect->Probability of individual nodes->The specific formula is as follows:
；
wherein the method comprises the steps ofIs->The number of elements;
calculation is given at the firstProbability of selecting conductivity values on individual tree nodes +.>The specific formula is as follows:
；
wherein,is a normal distribution of parameter disturbanceStandard deviation (S)>Is->Selecting a conductivity change value on each tree node;
the specific formula of the disturbance recommendation distribution in S6.3.2 is as follows:
；
the specific formula for generating the suggested distribution in S6.3.2 is as follows:
；
wherein,and->The probability of selecting the location where the new node is placed is expressed as follows:
；
wherein,the number of generatable node sets represented as the current model;
the calculation formula of the extinction proposal distribution in S6.3.2 is as follows:
；
wherein,is a set of nodes that may be deleted after the proposed birth,/for example>The number of sets of resolvable nodes expressed as candidate models.
Optionally, in S6:
selection of suggested model acceptance rates for existing nodesThe expression of (2) is:
；
selection of suggested model acceptance rates for generatable nodesThe expression of (2) is:
；
suggested model acceptance rate for selecting an evaluable nodeThe expression of (2) is:
；
wherein,number of set of resolvable nodes expressed as current model,/->The number of sets of nodes that can be generated, expressed as candidate models.
Optionally, in the step S7,，/>not less than 1000000%>、/>Taking natural number.
In order to realize a physical space tree structure, the invention improves the tree structure of a wavelet domain, and in the original tree structure, the node coefficient value is expressed as the wavelet coefficient without any physical meaning, so that we cannot accurately give the prior distribution condition of related coefficients, which leads the conductivity value obtained by inversion to be possibly beyond the actual condition, therefore, the invention changes the wavelet domain tree structure into the physical space tree structure, and the node parameter of the tree structure is the prior distribution with the electrical distribution characteristic. According to the invention, a likelihood function of the suggestion model is obtained according to the actually measured underground electromagnetic field response and the forward modeling response calculation of the suggestion model, and the acceptance rate of the suggestion model is calculated, so that the suggestion model is screened to obtain a sampling model which is closer to the actual underground conductivity model, and the conductivity parameters of the sampling model can be changed by optional nodes, so that inversion parameters can be flexibly adjusted, and the sampling can be conveniently carried out for a plurality of times; the invention solves the problem that the prior distribution is difficult to select by the nodes of the wavelet domain tree structure by adopting the physical space tree structure, achieves the effects of automating parameters and multi-scale subdivision of inversion areas, saves the search of a Markov chain on a low probability area, effectively improves the sampling efficiency, can be extended to Gao Weibei leaf inversion, and provides an automatic parameterized thought for a Gao Weibei leaf inversion method.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings that are required in the embodiments or the description of the prior art will be briefly described, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and other drawings may be obtained according to the structures shown in these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a design method of a cross-dimensional Bayesian sampler based on a physical space tree structure in an embodiment of the present invention;
FIG. 2 is a schematic diagram of a sampling flow in an embodiment of the invention;
FIG. 3 is a diagram of three different types of nodes included in a tree structure used in the present invention;
FIG. 4 shows a different arrangement of the same nodes in the embodiment of the present invention, wherein: fig. 4 (a) to 4 (e) show all tree structure arrangements when tree nodes (black boxes are shown) are 3;
FIG. 5 is a diagram illustrating a one-dimensional electrical structure corresponding to a binary tree in accordance with an embodiment of the present invention, wherein: fig. 5 (a) shows the conductivity distribution of the underground uniform half space corresponding to the binary tree when only the root node exists, fig. 5 (b) shows the underground conductivity distribution corresponding to the case where the first two layers of nodes exist in the tree structure, and fig. 5 (c) shows the underground conductivity distribution corresponding to the case where the first three layers of nodes exist in the tree structure;for a measured minimum subsurface depth,for a measured maximum subsurface depth;
FIG. 6 is a diagram of mapping functions used in the conversion from a tree structure to a physical model in the practice of the present invention, wherein: fig. 6 (a), 6 (b), 6 (c) show graphs of scale functions at different scale, fig. 6 (d), 6 (e), 6 (f) show graphs of decoupling functions at different scale (the first number of subscripts indicates scale, the second number indicates scale);
FIG. 7 is a prior verification diagram of tree nodes in an embodiment of the present invention, wherein: fig. 7 (a), fig. 7 (b), and fig. 7 (c) show the tree structure as a binary tree,Poise and parkThe loose coefficients are +.>A result graph at the time;
FIG. 8 is a graph of a one-dimensional edge probability density distribution (black line is a true model distribution) obtained by inversion in the practice of the present invention, wherein: the abscissa is the distribution of conductivity values, and the ordinate is the depth;
FIG. 9 is a graph of the distribution of RMS during sampling in an embodiment of the invention, wherein: the abscissa is the number of times of sampling, and the ordinate is the fitting difference of the data;
FIG. 10 is a graph of a sample node number distribution in an embodiment of the present invention, wherein: the abscissa is the node number, and the ordinate is the frequency;
fig. 11 is a result of sampling using a wavelet based tree structure, wherein: fig. 11 (a) shows 100 tens of thousands of times, fig. 11 (b) shows 500 tens of thousands of times, and fig. 11 (c) shows 1500 tens of thousands of times;
FIG. 12 is a graph showing the results of sampling using a physical space tree-based cross-dimensional Bayesian sampler in accordance with the present invention, in which: fig. 12 (a) shows the total number of samples as 100 ten thousand times, fig. 12 (b) shows the total number of samples as 500 ten thousand times, and fig. 12 (c) shows the result of the total number of samples as 1500 ten thousand times.
The achievement of the objects, functional features and advantages of the present invention will be further described with reference to the accompanying drawings, in conjunction with the embodiments.
Detailed Description
Embodiments of the invention are described in detail below with reference to the attached drawings, but the invention can be implemented in a number of different ways, which are defined and covered by the claims.
The invention provides a design method of a cross-dimensional Bayes sampler based on a physical space tree structure, and aims to solve the problems that prior distribution selection of nodes of a wavelet domain tree structure is difficult and sampling efficiency is low.
As shown in fig. 1, the design method of the cross-dimensional bayesian sampler based on the physical space tree structure comprises the following steps:
s1: measuring by using an electromagnetic measuring instrument to obtain the response of an underground electromagnetic field of a certain area, and determining the measured underground depth;
s2: constructing a tree structure initial model, constructing a Markov chain, and inputting data, the initial model and prior information as shown in fig. 2; calculating the fitting difference of the initial model, starting sampling,from 1 to a maximum number of samples;
the details are as follows:
calculating the maximum depth of the tree structure, wherein the formula is as follows:wherein->Representing the number of subdivisions of the subsurface depth at the maximum depth of the tree-like structure,/->Representing the maximum depth of the tree structure and constructing the maximum depth +.>A lower binary tree structure;
number of times of construction samplingTime tree structure initial model->Wherein->For the number of nodes in the initial model, +.>Depth level for tree structure, ++>， />Expressed as depth->Lower part (C)The number of arrangement of tree structures of all different shapes composed of the number of nodes, < >>Indicate->The conductivity parameters formed by the individual nodes are shown in fig. 5;
from the slaveWhen starting to build Markov chain for +.>Existing tree structure at subsampling +.>In other words, the nodes of the tree structure can be classified into existing nodes, generatable nodes, and deceased nodes;
s3: optional tree structureTree structure change is carried out on the class of nodes of the model (1) to form a new candidate model +.>The method comprises the steps of carrying out a first treatment on the surface of the If a generatable node is selected, a new candidate tree node is represented + ->The new candidate model isThe method comprises the steps of carrying out a first treatment on the surface of the If an erasable node is selected, a new candidate tree node is indicated +.>The new candidate model is +.>The method comprises the steps of carrying out a first treatment on the surface of the If one of the existing nodes is selected, it means that only the conductivity parameter of the selected node is changed +.>The new candidate model is +.>；
S4: new candidate models under the tree structure obtained in the step S3Performing inverse transformation to obtain a proposed model of the conductivity structure>And by forward algorithm ∈ ->Will suggest a model->Mapping to a data space to obtain forward response;
specifically, as shown in fig. 6, the inverse transformation in S4 specifically includes:
construction of a Scale functionTo control the currently existing nodes, the specific formula is as follows:
；
and for tree structure depth isLower->The individual nodes may define:
；
wherein,expressed as +.>Unit function on->Any number expressed as 0-1 is not specifically defined, ">Expressed as the level at which the node is located, +.>The position of the level of the node is shown;
constructing a decoupling functionThe relationship between the upper level node and the next level node is decoupled, namely the value of the parent node of the range where the child node is located is covered by the value of the child node, and the specific formula is as follows:
；
similarly, for tree structure depth isLower->The individual nodes define:
；
wherein,expressed as +.>Unit function on the table.
S5: calculating a suggestion model based on the electromagnetic field response in S1 and the forward response obtained in S4Likelihood functions of (2);
the method specifically comprises the following steps:
s5.1: suggested model obtained according to inverse transformation of tree-structured candidate modelLeast square method is carried out on the forward response of the model (A) and the electromagnetic field response obtained by measurement, and fitting difference +.>The specific formula is as follows:
；
wherein,representing the mapping of the model vector to the positive operator in the data space, < >>Is the data covariance matrix,>response data for the actual measured electromagnetic field;
s5.2: assuming that likelihood functions obey a standard normal distribution, likelihood functionsThe distribution of (2) is specifically expressed as follows:
；
wherein,for the number of observations.
S6: obtaining acceptance rate of suggestion model based on likelihood function calculation in S5Taking random number between 0 and 1, and judging the acceptance rate +.>Whether the random number is smaller than the random number, if yes, receiving a new candidate model to obtain the +.>The existing tree structure at sub-sampling is: />Otherwise, rejecting the new candidate model to obtain +.>The existing tree structure at sub-sampling is:；
the method specifically comprises the following steps:
s6.1: computing a priori distribution of a current existing tree structure given a binary tree structurePrior distribution of new candidate models +.>；
S6.2: calculate the firstExisting tree-like structure model at subsampling +.>Inverse transformation of the resulting conductivity model->Likelihood function of +.>First->New candidate tree-like structure model at subsampling +.>The proposed model obtained by inverse transformation->Likelihood function of +.>；
S6.3: calculation from existing tree structure modelTo a new candidate tree-structure model->Is a suggested distribution probability of (2)From the new candidate tree structure model +.>To the existing tree structure model->Is->；
S6.4: calculating the acceptance rate of the proposal model according to the prior distribution, likelihood function and proposal distribution calculated in S6.1 to S6.3The specific formula is as follows:
；
wherein,is Jacobian matrix->；
S6.5: taking a random number between 0 and 1, and judging the acceptance rateWhether or not the random number is smaller than the random number, if yes, accepting the advice model, i.e. +.>Otherwise reject the advice model, i.e. +.>。
S6.1 includes:
s6.1.1: calculating uniform distribution over a given range of node numbersThe specific formula is as follows:
；
wherein,is the maximum node number of the tree structure, +.>，/>Is the minimum node number of the tree structure, +.>Generally set to 1, indicating at least one node is owned;
s6.1.2: calculating the prior number of the arrangement under the given node numberNamely, the distribution of the permutation and combination number of the tree structure under the condition of a given node is as follows:
；
wherein,is->The number of nodes forms the number of tree-like structures of all different shapes,，/>is a standard binomial coefficient;
in one embodiment, a given number of nodesThe tree structure arrangement may be as shown in fig. 4.
S6.1.3: computing a priori distribution of parameters in each given nodeThe specific formula is as follows:
；
s6.1.4: the prior distribution under the binary tree structure is calculated respectively by adopting the same formulaAnd->Wherein a priori distribution->The specific formula of (2) is as follows:
。
s6.3 comprises:
s6.3.1: forming the current nodes on the current tree structure into a perturbable node groupThe generatable nodes form a generatable node group +.>The eroding nodes form eroding node group +.>；
As shown in FIG. 3, the resolvable node group represents all nodes of the child nodes which are not active on the current tree structure, the perturbable node group represents all nodes of the child nodes which are not active on the current tree structure and only one child node, and the generatable node group is a collection of hollow nodes in the tree structure.
S6.3.2: if the node selected in the S3 is the existing node, calculating disturbance suggestion distribution, if the node selected in the S3 is a generatable node, calculating generation suggestion distribution, and if the node selected in the S3 is an deceaseable node, calculating deceased suggestion distribution.
The calculation of the disturbance recommendation distribution in S6.3.2 includes:
calculating a givenSelect->Probability of individual nodes->The specific formula is as follows:
；
wherein the method comprises the steps ofIs->The number of elements;
calculation is given at the firstProbability of selecting conductivity values on individual tree nodes +.>The specific formula is as follows: />
；
Wherein,standard deviation of normal distribution for parameter perturbation, +.>First->Selecting a conductivity change value on each tree node;
the specific formula for calculating the disturbance recommendation distribution in S6.3.2 is as follows:
；
the specific formula for generating the suggested distribution in S6.3.2 is as follows:
；
wherein,and->The probability of selecting the location where the new node is placed is expressed as follows:
；
wherein,the number of generatable node sets represented as the current model;
the calculation formula of the extinction proposal distribution in S6.3.2 is as follows:
；
wherein,is a set of nodes that may be deleted after the proposed birth,/for example>The number of sets of resolvable nodes expressed as candidate models.
S5:
selection of suggested model acceptance rates for existing nodesThe expression of (2) is:
；
selection of suggested model acceptance rates for generatable nodesThe expression of (2) is:
；
suggested model acceptance rate for selecting an evaluable nodeThe expression of (2) is: />
；
Wherein,number of set of resolvable nodes expressed as current model,/->The number of sets of nodes that can be generated, expressed as candidate models.
S7: order theRepeating S3 to S6 and performing +/every time>Sampling model obtained by sampling once is stored until sampling times +.>And ending the time. Wherein (1)>，/>Not less than 1000000%>、/>Taking a natural number;
s8: and (3) carrying out statistical calculation on all the sampling models obtained in the step (S7), and obtaining the one-dimensional edge probability density distribution of the underground conductivity of the area, as shown in figure 8.
For such a tree-structured cross-dimensional sampler, it is critical whether the M-H acceptance criteria constructed above are correct. While a key check on the correctness of the acceptance criteria for a tree-structured cross-dimensional sampler is to test the number of nodes for the criteriaWhether posterior bias exists. One effective method is that the sampler operates a large number of steps, and regards likelihood function items as constant items, and finally compares the distribution condition of the sampled node number with the estimated distribution condition to judge whether the distribution accords with the expectations. For the likelihood function, a constant is maintained so that the perturbation step is not performed in the cross-dimensional sampler, the likelihood ratio in the generation step and the extinction step is kept constant at 1, and the probability of performing the two steps is maintained。
In the present embodiment, the number of nodes is100 ten thousand samples were tested, +.>For a maximum node value set by person, +.>. Due to->Is of poisson distribution, and the number of nodes can be obtained before sampling>Priori->The following is shown:
。
setting up10,20,30, respectively, the results of which are shown in fig. 7. The black bar graph in the graph represents the distribution of the number of nodes after sampling, and the black line represents the prior distribution to which the nodes should obey. From the figure it can be seen +.>=10，/>=20>The distribution obtained by sampling at=30 is quite consistent with the a priori known distribution.
In order to verify inversion effect of the method provided by the invention, thus proving that the required uncertainty information can be obtained, performing model inspection, designing an underground conductivity model shown in table 1, sampling by a designed sampler to obtain a posterior probability model, and describing in detail as follows:
table 1: conductivity and interface depth information
/>
The cross-dimension sampler designed by the invention is realized by utilizing Fortran language programming, and a personal notebook computer used for running a program is configured as follows: CPU-InterCore i7-8700, main frequency 3.2GHz, running memory 8.00GB. FIG. 8 is a graph of a one-dimensional edge probability density distribution of a distribution of an underground medium sampled by a sampler according to the present invention, wherein the shade of the edge probability density distribution indicates the conductivity probability of a certain region at that depth, and the darker the shade indicates the higher the probability. From the results obtained by inversion, the 5-layer model is better recovered. The first, second, third and last layers have better edge probability density constraint and small uncertainty, so that the sampling efficiency of the sampler for the layers is high, the uncertainty range of the fourth layer is large, the fact that the magnetotelluric method is insensitive to high resistance and the expressive force of the magnetotelluric method is supposed to have physical characteristics is shown, but from the probability of the magnetotelluric method, the high probability area is better corresponding to the real situation of the model, and the method can well extract the characteristic response of the high resistance layer in the data. And the sampling by the samplers shown in fig. 9 and 10 has good mixing and constraint on the number of parameters.
When inverting a conductivity model of three subsurface layers, when sampling using a wavelet based tree structure, it can be seen that the resulting distribution in FIG. 11 does not match the true model distribution until 500 tens of thousands of times after sampling. When the cross-dimensional Bayesian sampler based on the physical space tree structure is used for sampling, the distribution of the real model can be searched by only carrying out one million times of sampling in the figure 12. Therefore, the sampling efficiency is remarkably improved.
The foregoing description is only of the preferred embodiments of the present invention and is not intended to limit the scope of the invention, and all equivalent structural changes made by the description of the present invention and the accompanying drawings or direct/indirect application in other related technical fields are included in the scope of the invention.

Claims (9)

1. A design method of a cross-dimensional Bayesian sampler based on a physical space tree structure is characterized by comprising the following steps:
s1: measuring by using an electromagnetic measuring instrument to obtain the response of an underground electromagnetic field of a certain area, and determining the measured underground depth;
s2: constructing a tree structure initial model, and constructing a Markov chain, wherein the method specifically comprises the following steps of:
calculating the maximum depth of the tree structure, wherein the formula is as follows:wherein->Representing the maximum depth of the tree-like structure, +.>Representing the number of subdivisions of the subsurface depth at the maximum depth of the tree structure and constructing a maximum depth +.>A lower binary tree structure;
number of times of construction samplingTime tree structure initial model->Wherein->For the number of nodes in the initial model, +.>Depth level for tree structure, ++>， />Expressed as depth->Lower->The number of arrangement of tree structures of all different shapes composed of the number of nodes, < >>Indicate->Conductivity parameters formed by the individual nodes;
from the slaveWhen starting to build Markov chain for +.>Existing tree structure at subsampling +.>In other words, the nodes of the tree structure are classified into existing nodes, generatable nodes and deceased nodes;
s3: optional tree structureTree structure change is carried out on the class of nodes of the model (1) to form a new candidate model +.>The method comprises the steps of carrying out a first treatment on the surface of the If a generatable node is selected, a new candidate tree node is represented + ->The new candidate model is +.>The method comprises the steps of carrying out a first treatment on the surface of the If an erasable node is selected, a new candidate tree node is indicated +.>The new candidate model isThe method comprises the steps of carrying out a first treatment on the surface of the Indicating that only the conductivity parameter of the selected node is changed if one of the existing nodes is selectedNew candidate modeIs->；
S4: new candidate models under the tree structure obtained in the step S3Performing inverse transformation to obtain a suggested model of the conductivity structureAnd by forward algorithm ∈ ->Will suggest a model->Mapping to a data space to obtain forward response;
s5: calculating a suggestion model based on the electromagnetic field response in S1 and the forward response obtained in S4Likelihood functions of (2);
s6: obtaining acceptance rate of suggestion model based on likelihood function calculation in S5Judging the acceptance rate->Whether the random number is smaller than the random number, if yes, receiving a new candidate model to obtain the +.>The existing tree structure at sub-sampling is: />Otherwise, rejecting the new candidate model to obtain +.>The existing tree structure at sub-sampling is: />；
S7: order theRepeating S3 to S6 and performing +/every time>Sampling model obtained by sampling once is stored until sampling times +.>Ending the time;
s8: and (3) carrying out statistical calculation on all the sampling models obtained in the step (S7) to obtain the one-dimensional edge probability density distribution of the underground conductivity of the area.
2. The method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 1, wherein the inverse transformation in S4 specifically comprises:
construction of a Scale functionTo control the currently existing nodes, the specific formula is as follows:
；
and for tree structure depth isLower->The individual nodes are defined as:
；
wherein:expressed as +.>Unit function on; />Any number represented as 0 to 1 has no specific meaning; />Expressed as the level at which the node is located; />The position of the level of the node is shown;
constructing a decoupling functionThe relationship between the upper level node and the next level node is decoupled, namely the value of the parent node of the range where the child node is located is covered by the value of the child node, and the specific formula is as follows:
；
similarly, for tree structure depth isLower->The individual nodes are defined as:
；
wherein,expressed as +.>Unit function on the table.
3. The method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 1, wherein S5 comprises:
s5.1: suggested model obtained according to inverse transformation of tree-structured candidate modelLeast square method is carried out on the forward response of the model (A) and the electromagnetic field response obtained by measurement, and fitting difference +.>The specific formula is as follows:
；
wherein,representing the mapping of the model vector to the positive operator in the data space, < >>Is the data covariance matrix,>response data for the actual measured electromagnetic field;
s5.2: assuming that likelihood functions obey a standard normal distribution, likelihood functionsIs provided with (1)The volume formula is as follows:
；
wherein,for the number of observations.
4. A method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 3 and wherein S6 comprises:
s6.1: computing a priori distribution of a current existing tree structure given a binary tree structurePrior distribution of new candidate models +.>；
S6.2: calculate the firstExisting tree-like structure model at subsampling +.>Inverse transformation of the resulting conductivity model->Likelihood function of +.>First->New candidate tree-like structure model at subsampling +.>Through inverse transformation to obtainSuggested model->Likelihood function of +.>；
S6.3: calculation from existing tree structure modelTo a new candidate tree-structure model->Is a suggested distribution probability of (2)From the new candidate tree structure model +.>To the existing tree structure model->Is->；
S6.4: calculating the acceptance rate of the proposal model according to the prior distribution, likelihood function and proposal distribution calculated in S6.1 to S6.3The specific formula is as follows:
；
wherein,is Jacobian matrix->；
S6.5: taking a random number between 0 and 1, and judging the acceptance rateWhether or not the number is smaller than the random number, if so, accepting the proposal model, namelyOtherwise reject the advice model, i.e. +.>。
5. The method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 4 wherein the step S6.1 comprises:
s6.1.1: calculating uniform distribution over a given range of node numbersThe specific formula is as follows:
；
wherein,is the maximum node number of the tree structure, +.>；/>Is the minimum node number of the tree structure, +.>Generally set to 1, indicating at least one node is owned;
S6.1.2: calculating the prior number of the arrangement under the given node numberNamely, the distribution of the permutation and combination number of the tree structure under the condition of a given node is as follows:
；
wherein,is->The number of nodes forms the number of tree-like structures of all different shapes,，/>is a standard binomial coefficient;
s6.1.3: computing a priori distribution of parameters in each given nodeThe specific formula is as follows:
；
s6.1.4: the prior distribution under the binary tree structure is calculated respectively by adopting the same formulaAnd->Wherein a priori distribution->The specific formula of (2) is as follows:
。
6. the method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 5 wherein S6.3 comprises:
s6.3.1: forming the current nodes on the current tree structure into a perturbable node groupThe generatable nodes form a generatable node group +.>The eroding nodes form eroding node group +.>；
S6.3.2: if the node selected in the S3 is the existing node, calculating disturbance suggestion distribution, if the node selected in the S3 is a generatable node, calculating generation suggestion distribution, and if the node selected in the S3 is an deceaseable node, calculating deceased suggestion distribution.
7. The method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 6, wherein the calculating of the disturbance recommendation distribution in S6.3.2 comprises:
calculating a givenSelect->Probability of individual nodes->The specific formula is as follows:
；
wherein the method comprises the steps ofIs->The number of elements;
calculation is given at the firstProbability of selecting conductivity values on individual tree nodes +.>The specific formula is as follows:
；
wherein,standard deviation of normal distribution for parameter perturbation, +.>Is->Conductivity change values at individual tree nodes;
the specific formula of the disturbance recommendation distribution in S6.3.2 is as follows:
；
the specific formula for generating the suggested distribution in S6.3.2 is as follows:
；
wherein,and->The probability of selecting the location where the new node is placed is expressed as follows:
；
wherein,the number of generatable node sets represented as the current model;
the calculation formula of the extinction proposal distribution in S6.3.2 is as follows:
；
wherein,is a set of nodes that may be deleted after the proposed birth,/for example>The number of sets of resolvable nodes expressed as candidate models.
8. The method for designing a cross-dimensional bayesian sampler based on a physical space tree structure according to claim 7, wherein in S6:
selection of suggested model acceptance rates for existing nodesThe expression of (2) is:
；
selection of suggested model acceptance rates for generatable nodesThe expression of (2) is:
；
suggested model acceptance rate for selecting an evaluable nodeThe expression of (2) is:
；
wherein,number of set of resolvable nodes expressed as current model,/->The number of sets of nodes that can be generated, expressed as candidate models.
9. A method for designing a cross-dimensional Bayesian sampler based on a physical space tree structure as claimed in any of claims 1-8, wherein in S7,，/>not less than 1000000%>、/>Taking natural number.
CN202410078165.4A 2024-01-19 2024-01-19 Design method of cross-dimension Bayesian sampler based on physical space tree structure Active CN117574790B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202410078165.4A CN117574790B (en) 2024-01-19 2024-01-19 Design method of cross-dimension Bayesian sampler based on physical space tree structure

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202410078165.4A CN117574790B (en) 2024-01-19 2024-01-19 Design method of cross-dimension Bayesian sampler based on physical space tree structure

Publications (2)

Publication Number Publication Date
CN117574790A true CN117574790A (en) 2024-02-20
CN117574790B CN117574790B (en) 2024-03-26

Family

ID=89864865

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202410078165.4A Active CN117574790B (en) 2024-01-19 2024-01-19 Design method of cross-dimension Bayesian sampler based on physical space tree structure

Country Status (1)

Country Link
CN (1) CN117574790B (en)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010151354A1 (en) * 2009-06-26 2010-12-29 Exxonmobil Upstream Research Company Constructing resistivity models from stochastic inversion
US20190011583A1 (en) * 2017-07-06 2019-01-10 Chevron U.S.A. Inc. System and method for full waveform inversion of seismic data
CN113177330A (en) * 2021-05-27 2021-07-27 吉林大学 Transient electromagnetic rapid statistical inversion method
CN113553773A (en) * 2021-08-16 2021-10-26 吉林大学 Ground-air electromagnetic data inversion method based on Bayesian framework combined with neural network
CN114065511A (en) * 2021-11-15 2022-02-18 中南大学 Magnetotelluric two-dimensional forward modeling numerical simulation method, device, equipment and medium under undulating terrain
CN114896564A (en) * 2022-05-23 2022-08-12 武汉市市政建设集团有限公司 Transient electromagnetic two-dimensional Bayesian inversion method adopting self-adaptive Thiessen polygon parameterization
CN116992668A (en) * 2023-08-03 2023-11-03 中国科学院地质与地球物理研究所 Transient electromagnetic self-adaptive transverse constraint inversion method

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010151354A1 (en) * 2009-06-26 2010-12-29 Exxonmobil Upstream Research Company Constructing resistivity models from stochastic inversion
US20190011583A1 (en) * 2017-07-06 2019-01-10 Chevron U.S.A. Inc. System and method for full waveform inversion of seismic data
CN113177330A (en) * 2021-05-27 2021-07-27 吉林大学 Transient electromagnetic rapid statistical inversion method
CN113553773A (en) * 2021-08-16 2021-10-26 吉林大学 Ground-air electromagnetic data inversion method based on Bayesian framework combined with neural network
CN114065511A (en) * 2021-11-15 2022-02-18 中南大学 Magnetotelluric two-dimensional forward modeling numerical simulation method, device, equipment and medium under undulating terrain
CN114896564A (en) * 2022-05-23 2022-08-12 武汉市市政建设集团有限公司 Transient electromagnetic two-dimensional Bayesian inversion method adopting self-adaptive Thiessen polygon parameterization
CN116992668A (en) * 2023-08-03 2023-11-03 中国科学院地质与地球物理研究所 Transient electromagnetic self-adaptive transverse constraint inversion method

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
田圣琦: "基于小波树状结构的二维大地电磁跨维贝叶斯反演", 中国优秀硕士学位论文全文数据库 基础科学辑, no. 01, 15 January 2024 (2024-01-15), pages 011 - 301 *

Also Published As

Publication number Publication date
CN117574790B (en) 2024-03-26

Similar Documents

Publication Publication Date Title
Batista et al. On the evaluation of soil erosion models: Are we doing enough?
Al-Bulushi et al. Artificial neural networks workflow and its application in the petroleum industry
Ma et al. Integration of artificial intelligence and production data analysis for shale heterogeneity characterization in steam-assisted gravity-drainage reservoirs
Demyanov et al. Neural network residual kriging application for climatic data
de Fondeville et al. Functional peaks-over-threshold analysis
Li et al. Novel multiple resolutions design of experiment/response surface methodology for uncertainty analysis of reservoir simulation forecasts
CN113687433B (en) Bi-LSTM-based magnetotelluric signal denoising method and system
Loquin et al. A fuzzy interval analysis approach to kriging with ill-known variogram and data
Jardani et al. Use of convolutional neural networks with encoder-decoder structure for predicting the inverse operator in hydraulic tomography
Le Goc et al. An inverse problem methodology to identify flow channels in fractured media using synthetic steady-state head and geometrical data
Brantson et al. Development of machine learning predictive models for history matching tight gas carbonate reservoir production profiles
Rotondi et al. Hydrocarbon production forecast and uncertainty quantification: A field application
Chen et al. NMR-data-driven prediction of matrix permeability in sandstone aquifers
Roudini et al. Proxy-based Bayesian inversion of strain tensor data measured during well tests
Guo et al. Groundwater depth forecasting using configurational entropy spectral analyses with the optimal input
Miao et al. BayLUP: A Bayesian framework for conditional random field simulation of the liquefaction-induced settlement considering statistical uncertainty and model error
Yanshu et al. A three-dimensional model of deep-water turbidity channel in Plutonio oilfield, Angola: From training image generation, optimization to multi-point geostatistical modelling
Yu et al. Gated recurrent unit neural network (GRU) based on quantile regression (QR) predicts reservoir parameters through well logging data
CN117574790B (en) Design method of cross-dimension Bayesian sampler based on physical space tree structure
Gallop Facies probability from mixture distributions with non-stationary impedance errors
Morton et al. Grid-based inversion methods for spatial feature identification and parameter estimation from pressure transient tests
EP3929631B1 (en) Method and system for analyzing a reservoir grid of a reservoir geological formation based on 4d seismic images
Heydari Gholanlo et al. Estimation of cementation factor in carbonate reservoir by using genetic fuzzy inference system
Ye et al. Rapid and Consistent Identification of Stratigraphic Boundaries and Stacking Patterns in Well Logs-An Automated Process Utilizing Wavelet Transforms and Beta Distributions
Thiam et al. Reservoir interwell connectivity estimation from small datasets using a probabilistic data driven approach and uncertainty quantification

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant