WO2011159255A2 - Analyse de données à grand nombre de dimensions - Google Patents

Analyse de données à grand nombre de dimensions Download PDF

Info

Publication number
WO2011159255A2
WO2011159255A2 PCT/SG2011/000207 SG2011000207W WO2011159255A2 WO 2011159255 A2 WO2011159255 A2 WO 2011159255A2 SG 2011000207 W SG2011000207 W SG 2011000207W WO 2011159255 A2 WO2011159255 A2 WO 2011159255A2
Authority
WO
WIPO (PCT)
Prior art keywords
model parameter
parameter set
reduced base
sampling
base
Prior art date
Application number
PCT/SG2011/000207
Other languages
English (en)
Other versions
WO2011159255A3 (fr
Inventor
Juan Luis FERNÁNDEZ MARTINEZ
Original Assignee
Blue Prism Technologies Pte Ltd
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
Priority claimed from US13/041,423 external-priority patent/US8688616B2/en
Application filed by Blue Prism Technologies Pte Ltd filed Critical Blue Prism Technologies Pte Ltd
Publication of WO2011159255A2 publication Critical patent/WO2011159255A2/fr
Publication of WO2011159255A3 publication Critical patent/WO2011159255A3/fr

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V20/00—Geomodelling in general
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/32—Transforming one recording into another or one representation into another
    • G01V1/325—Transforming one representation into another
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/40—Transforming data representation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface

Definitions

  • the present disclosure relates to data analysis and, more specifically, to systems and methods for analyzing data in high-dimensional space.
  • Inverse problems are often encountered in many fields of technology, including engineering, science and mathematics. Solving an inverse problem entails the determination of certain model parameters from a set of observed data or measurements.
  • the mapping of data to model parameters is a result of interactions within a physical system, such as the Earth, atmosphere, gravity, and so forth.
  • geologic model parameters e.g., conductivity, density, magnetic permeability, porosity, seismic velocity, etc.
  • m (mi, n3 ⁇ 4 e M c R" denotes a model parameter set that belongs to a set or family of admissible model parameters M
  • d e R represents the observed data
  • F(m) (fi (m), f 2 (m), f s (m)) represents the forward model that predicts the observed data.
  • Uncertainty may be defined as the difference between the one true value that describes a physical quantity at a specific point in time and space and the value reported as a result of a measurement. Estimation of uncertainty involves finding the family M of equivalent model parameter sets m that are consistent with the prior knowledge and fit the observed data d e R s within a prescribed tolerance ( ⁇ ), as follows:
  • Quantifying uncertainly is a key aspect in risk management, business analysis, probabilistic forecasting and many other business processes.
  • Model-based applications that incorporate uncertainty evaluation capabilities can provide invaluable guidance in business decision-making, such as whether or not to acquire more data to reduce uncertainty or to proceed with the current path, or whether the potential reward that can be achieved in developing a set of assets is outweighed by the degree of risk and cost involved. For example, in the petrochemical field, the success of finding new oil and gas reserves can be significantly improved by evaluating the uncertainty of developing new leads.
  • Bayesian network-based frameworks have previously been used to estimate uncertainty.
  • such frameworks are very inefficient, especially where the number of parameters is very large and/or the forward evaluations are very costly to compute.
  • Sampling within a Bayesian framework incurs very high computational costs because it is performed in parts of the model space with very small likelihood of being consistent with the observed data.
  • these conventional methods are naturally limited to small parameterizations (i.e. low number of parameters) and fast forward solvers.
  • observed data and at least one input model parameter set is received.
  • the input model parameter set serves as a solution candidate of a predefined problem (e.g., inverse or optimization problem) and is related to the observed data via a model.
  • a reduced base with lower dimensionality is determined based on the input model parameter set.
  • the reduced base is associated with a set of coefficients, which represents the coordinates of any model parameter set in the reduced base. Sampling is performed within the reduced base to generate an output model parameter set in the reduced base that is compatible with the input model parameter set and fits the observed data, via the model, within a predetermined threshold.
  • FIG. 1 shows an exemplary system
  • FIG. 2 shows an exemplary method
  • FIG. 3a shows an exemplary method of constructing an orthogonal reduced base
  • FIG. 3b shows images representing exemplary reduced base terms with increasing spatial frequencies
  • FIG. 4 shows an exemplary application of singular value decomposition on a model
  • FIG. 5 shows an exemplary application of discrete cosine transform on a model
  • FIG. 6 shows an exemplary application of discrete wavelet transform on a model
  • FIG. 7 shows an exemplary application of discrete Chebyshev transform on a model
  • FIG. 8 shows an exemplary application of independent component analysis on a model
  • FIG. 9 shows another exemplary method for computing a spatial indicator reduced base
  • FIG. 10 shows an exemplary application of spatial indicator-based method on a model
  • FIG. 1 1 shows an exemplary bounded constrained sampling method
  • FIG. 12 shows an exemplary telescopic sampling method
  • FIG. 13 shows an exemplary plot of the valleys of equivalent model parameters
  • FIG. 14a shows an exemplary method of computing a posterior reduced base after optimization in a prior reduced base
  • FIG. 14b shows an exemplary method of stochastic sampling in the posterior reduced base
  • FIG. 15 shows an exemplary method of sampling while optimizing in the reduced base
  • FIG. 16 shows a 2-D plot generated as a result of applying an exemplary method of sampling while optimizing for a valley shape.
  • the present framework includes building a lower- dimensional reduced base over which sampling for equivalent models may be efficiently performed.
  • the reduced base is associated with a smaller set of coefficients than the input model parameter set(s).
  • the coefficients of the reduced base represent the coordinates of any model parameter set in the reduced base.
  • the reduced base may be constructed using orthogonal transformations, other model reduction techniques or any combination thereof.
  • the model reduction techniques include, but are not limited to, Singular Value Decomposition, Discrete Cosine Transform, Discrete Wavelet Transform, Discrete Chebyshev Transform, and Independent Component Analysis, Indicator-base (or set level base) technique, or Principal Component Analysis (PCA) using unconditional prior scenarios.
  • sampling may be performed within the reduced base to generate one or more output model parameter sets that fit the observed data within a predetermined error tolerance and are compatible with the input model parameter set(s).
  • Sampling techniques include, for example, geometric sampling (e.g., hyperprism sampling), oriented sampling, sampling while optimizing in the reduced base, stochastic sampling, or a combination thereof.
  • the present framework also provides the flexibility of using other types of model parameterizations, such as logarithmic sampling, to better condition the sampling procedure.
  • the orthogonal nature of the reduced base also allows high frequencies to be added in the sampling as needed, such as in telescopic sampling.
  • the resulting output model parameter set may then be analyzed and used to construct various kinds of measures of interest, such as probability maps or uncertainty measures that are useful for risk analysis or other types of assessments.
  • the present framework casts the problem of estimating uncertainty in an optimization framework that makes use of the natural geometric and algebraic interpretation of the uncertainty problem.
  • the present framework is based on the assertion that equivalent model parameter sets are located along flat elongated valleys of the misfit function landscape.
  • a misfit function refers to any objective measure of the difference between predicted data and observed (or measured) data.
  • the misfit function may also exhibit different local optima located in the different basins of the misfit function landscape. Such multimodality feature is prevalent in many inverse problems.
  • equivalent model parameter sets are locally associated with a hyperquadric of equivalence with axes that are oriented along the singular vectors of the local Jacobian operator (J) and dimensions that are proportional to the inverse of the eigenvalues.
  • J Jacobian operator
  • £ tol 2 ,— (7)
  • Am ⁇ denotes the vector (m -m 0 ) referred to as the base V
  • Ad 3 ⁇ 4 is the error data prediction F (m 0 ) - d referred to as the base U.
  • Equation (7) may be simplified as follows:
  • V right singular vectors of the Jacobian
  • the axes are proportional to 1 / in each
  • the data is insensitive to model variations if the inverse problem is linear. Therefore, it can be deduced that the low misfit equivalent region is defined by a valley shape around the components that belong to the null-space of JF mo , with the valley having an infinite length. This means that along these directions, the data is completely insensitive to model variations if the inverse problem was linear.
  • the cost function landscape is locally either a hyper ellipsoid or an elliptical cylinder.
  • the nonlinearity of the forward operator causes the nonlinear equivalent regions to have a croissant or tube-like shape with either a finite or infinite length.
  • Nonlinear equivalent regions of infinite length embodied as meandering tubes in the model space, may be related to finite resolution issues of the forward model. For linear inverse problems, these valleys are straight and are related to the null space of the forward linear operator.
  • a scalar optimization problem may be formulated as: min E (x) , xeR" ' wherein E (x) : R" - R is an w-dimensional scalar field.
  • Solving a scalar optimization problem includes finding the family of model parameters x that fulfils the condition E (x) ⁇ tol . Assuming a model parameter set x 0 : E (x 0 ) ⁇ tol. , for any other model parameter set x in the neighbourhood of x 0 , the Taylor expansion of the scalar misfit function E(x) is as follows:
  • the matrix HE ⁇ j is symmetric, and therefore, exhibits orthogonal diagonalization and can be decomposed into the following form:
  • V is a matrix composed of the eigenvectors of HE, >
  • D is the diagonal matrix constructed from the corresponding eigenvalues
  • V T is the inverse matrix of V
  • ⁇ (r r) contains the r non-null eigenvalues of D.
  • the last n-r eigenvectors contained in V span the null space of the Hessian in x 0 . Accordingly, the equation of the hyperquadric of equivalence is as follows:
  • the local hyperquadric of equivalence can change its type from ellipsoids (positive eigenvalues), hyperboloids (negative eigenvalues), or elliptical or hyperbolic cylinders (null eigenvalues).
  • the Hessian may have a semi-definite positive character and the local hyperquadric is either an ellipsoid or an elliptical cylinder.
  • the centre of the hyperquadric coincides with the Gauss-Newton solution of the optimization problem when the input model is located in x 0 .
  • embodiments of the present invention can be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof.
  • the present technology can be
  • the application program can be uploaded to, and executed by, a machine comprising any suitable architecture.
  • the system and method of the present disclosure may be implemented in the form of a software application running on a computer system, for example, a laptop, personal computer (PC), workstation, client device, mini-computer, storage system, handheld computer, server, mainframe computer, dedicated digital appliance, and so forth.
  • the software application may be stored on a non-transitory media locally accessible by the computer system and accessible via a hard wired or wireless connection to a network, for example, a local area network, or the Internet.
  • FIG. 1 shows an exemplary system 100 operable to implement a method and system of the present disclosure.
  • the system 100 may include a computer system 102, which may be a server, a network of servers or resources, a desktop personal computer, a portable laptop computer, another portable device, a mini-computer, a mainframe computer, a storage system, a dedicated digital appliance, or another device having a storage sub-system configured to store a collection of digital data items.
  • a computer system 102 which may be a server, a network of servers or resources, a desktop personal computer, a portable laptop computer, another portable device, a mini-computer, a mainframe computer, a storage system, a dedicated digital appliance, or another device having a storage sub-system configured to store a collection of digital data items.
  • the computer system 102 includes a processor 101, a network controller 103, non-transitory computer-readable media 104, an input/output interface 105, one or more input devices 109 (e.g., keyboard, mouse, touch screen, tablet, etc), and one or more output devices 1 1 1 (e.g., monitor or display).
  • input devices 109 e.g., keyboard, mouse, touch screen, tablet, etc
  • output devices 1 1 1 e.g., monitor or display
  • Various other peripheral devices such as an additional data storage device and a printing device, can also be connected to the computer system 102.
  • Non-transitory computer-readable media 104 can include random access memory (RAM), read only memory (ROM), magnetic floppy disk, disk drive, tape drive, flash memory, etc., or a combination thereof.
  • the present invention may be implemented as a problem solver 120 and an uncertainly assessor 122 that includes computer-readable program code tangibly embodied in the non-transitory computer-readable media 104 and executed by the processor 101.
  • the computer system 102 is a general purpose computer system that becomes a specific purpose computer system when executing the routine of the present invention.
  • the computer-readable program code is not intended to be limited to any particular programming language and implementation thereof. It will be appreciated that a variety of programming languages and coding thereof may be used to implement the teachings of the disclosure contained herein.
  • the computer system 102 may also include an operating system and micro instruction code. The various processes and functions described herein can either be part of the micro instruction code or part of the application program or routine (or combination thereof) which is executed by the processor 101 via the operating system.
  • the computer system 102 may be connected, via network 132 using a communication link 130, to one or more client workstations 134 and data source 145.
  • the communication link may be a telephone, a wireless network link, a wired network link, or a cable network.
  • the network 132 may be wired or wireless, and may include a local area network (LAN), a wide area network (WAN), or a combination thereof.
  • the client workstation 134 may include a computer and appropriate peripherals, such as a keyboard and display, and can be operated in conjunction with the entire system 100.
  • the client workstation 134 may further include a memory to store information (e.g., input model parameter sets, prior scenarios, data, forward model definitions, etc.) that defines the inverse or optimization problem to be solved by the computer system 102.
  • the client workstation 134 may further include a user interface to allow a user (e.g., customer or client) to interactively manipulate the retrieved and/or processed data. For example, the user may specify the problem to be solved and provide instructions to the computer system 102.
  • the client workstation 134 may be communicatively coupled, via the network 132, with the data source 145 so that the data stored or collected by the data source 145 can be retrieved by the client workstation 134.
  • the data source 145 may be incorporated in the client workstation 134 or the computer system 102.
  • the data source 145 is a memory or other program storage device that stores digital observed data that can be used by the computer system 102 to predict model parameter sets.
  • the data source 145 may also be a data acquisition system that is adapted to collect the digital observed data by quantifying observations of a physical system. For example, such digital observed data may be geophysical data,
  • meteorological data or biomedical data is not limited to data observed in a physical system.
  • Non-physical systems such as financial markets or business enterprises, may also be observed.
  • the observed data may include historical stock prices, commodity prices, currency values or other types of financial data collected for analyzing investment portfolios and predicting future performance.
  • the observed data may be in any digital format, such as one-dimensional data or two- or three-dimensional images.
  • the data source 145 is a data acquisition system that includes a signal sensor, detector, receiver, source, scanner or any other suitable device that is operable to digitize observations (e.g., images) of a physical system (e.g., Earth, atmosphere, biological body, etc.).
  • the data acquisition system 145 may be adapted to collect geophysical data such as seismic reflection data that can be used to compute quantitative parameters (e.g., velocity structure, elastic parameters, etc.) of rock terrains via seismic inversion techniques.
  • seismic reflection techniques are particularly useful in hydrocarbon exploration, where subsurface distribution of terrains and its structure are mapped to delineate and assess the likelihood of potential hydrocarbon accumulations or ore deposits in oil and gas exploration.
  • Uncertainly in such data may exist due to, for example, partial coverage of the survey or measurement.
  • geophysical surveys may be incomplete because they are typically performed from the surface of the earth.
  • Uncertainty may also be caused by noise introduced during data acquisition, such as instrumental noise, noise caused by error in the location of sources and receivers, and so forth.
  • noise introduced during data acquisition such as instrumental noise, noise caused by error in the location of sources and receivers, and so forth.
  • biological data such as electrophysiological data or radiological images (e.g., magnetic resonance (MR) or computed tomography (CT) images) may be collected to detect abnormal medical conditions (e.g., tumor or malignant growths).
  • MR magnetic resonance
  • CT computed tomography
  • weather forecasting may be performed based on meteorological data indicative of precipitation, temperature, humidity values, cloud coverage, air quality, contamination dispersion, etc.
  • historical production and pressure data is collected to predict the flow of fluids (e.g., oil, water, gas) through porous media.
  • Relevant model parameters include, for example, permeability, porosity and saturation.
  • Other types of applications include biomolecular structural predictions, drug design, structural optimization, material design and nanotechnology, semiconductor design, chemo-metrics, and so forth.
  • the observed data (e.g., acquired or measured data) is transmitted from the data source 145 (or the client workstation 134) to the computer system 102 for processing.
  • Other types of information may also be sent to the computer system 102 for processing.
  • the client workstation 134 may provide an input model parameter set that is determined to be a possible solution candidate of a predefined problem (e.g., inverse or optimization problem).
  • the client workstation 134 may provide prior information (or prior scenarios) from which the computer system 102 may extract the input model parameter set.
  • the client workstation 134 may provide a mathematical model (e.g., forward model) that maps the input model parameter set to the observed data.
  • the problem solver 120 and uncertainty assessor 122 are executed by the processor 101 to process the observed data (e.g., measured or acquired data) and any input information provided.
  • the problem solver 120 generates a collection of output model parameter sets in a reduced base.
  • the output model parameter sets are compatible with the input model parameter set and fits the observed data within a predetermined error tolerance.
  • the uncertainty assessor 122 then processes the output model parameter sets to construct one or more uncertainty measures (e.g., probability map) that are useful for risk analysis or other types of uncertainty-related assessments.
  • the output model parameter set and/or results of the uncertainty assessment may be transmitted to the client workstation 134 for further processing or display.
  • FIG. 2 shows an exemplary method 200 of solving a predefined problem.
  • the method 200 may be implemented by, for example, the computer system 102 described with reference to FIG. 1.
  • the predefined problem may be an inverse problem, or a more general optimization problem, as discussed previously.
  • An inverse problem may be treated like an optimization problem and/or sampling problem.
  • the problem solver 120 receives observed data (e.g., measured data) from, for example, the data source 145.
  • the data source 145 may be located externally or incorporated in the client workstation 134 or computer system 102.
  • at least one input model parameter set may also be received.
  • the input model parameter set serves as a candidate solution and is related to the observed data via a model.
  • the model is a description of a system using mathematical language.
  • the model may be a forward model which predicts the observed data from parameters in the model space (or parameter space).
  • Mathematical models can take many forms, including but not limited to, dynamical systems, statistical models, a set of partial (or ordinary) differential equations, a set of algebraic equations, game theoretic models, etc. It is also understood that a model space, as used herein, refers to a multidimensional region describing the set of all possible values that the model parameters can take. Each point in the model space represents a model parameter set, which is a candidate or potential solution to the predefined problem. Each model parameter set includes a given number ( «) of unit values or model parameters (e.g., pixels of an image).
  • the input model parameter set comprises a unique posterior model parameter set mo that is located in a valley of multiple solutions.
  • the unique posterior model parameter set mo may be derived by using local optimization methods that converge to one of these multiple solutions.
  • an ensemble of input model parameter sets (mi, ⁇ 2, 1113, 7) may be derived from prior information (or prior scenarios). This is achieved by, for example, conditional and/or non-conditional simulation techniques, stochastic modeling or any other suitable techniques.
  • Conditional simulation techniques include, but are not limited to, multipoint geostatistical methods, turning bands, sequential simulation, simulated annealing, probability field simulation, matrix decomposition methods, and so forth.
  • the problem solver 120 computes a reduced base using the input model parameter set(s). This is performed in order to reduce the dimensionality of the correlated search space and create a reduced base associated with coefficients that are consistent with the structure of the input model parameter set and observed data. These coefficients represent the coordinates of any model parameter set located in the reduced base, and are much fewer than the number of model parameters in the input model parameter set, thereby rendering the subsequent search in the reduced base more efficient.
  • I is a space of rectangular matrices M(s,n) .
  • the canonical base associated with the input model parameter set has x n terms, but this is not very informative in guiding the search in the low misfit regions of the model space.
  • lower dimensional bases i.e. reduced bases associated with subspaces of the model parameter space I may be constructed.
  • Dimensionality of the reduced base may be established by, for example, assessing the percentage of model compression achieved. More particularly, the process of reducing the model base may be described as a process of finding a set of independent vectors ⁇ , , v 2 , ... , v ? ⁇ which are consistent with our prior knowledge.
  • the coordinates of a model parameter set m in the reduced base may be represented as follows:
  • ⁇ is a fixed vector in R" or trend (which can eventually be the null-vector) and (a, a q ) are coefficients of the model parameter set m in the reduced base.
  • the dimensionality of the reduced base is much lower, since the number q of coefficients a in the reduced base is much lower than the number n of the model parameters in the input model parameter set.
  • the reduced base is a local base situated along the valley of the misfit landscape that represents model parameters in an optimal way.
  • the reduced base may also be constructed from a set of prior scenarios (or input model parameter sets) using, for example, principal component analysis (PCA).
  • PCA principal component analysis
  • the reduced base is constructed based on the input model parameter set(s).
  • the reduced base may be a prior base or a posterior base.
  • the input model parameter set may include a single posterior model parameter set m 0 or an ensemble of model parameter sets derived from multiple prior scenarios.
  • a posterior reduced base may be built.
  • a prior reduced base may be constructed by orthogonalizing these prior scenarios using, for example, PCA or kernel PCA.
  • the prior or posterior reduced base may be built using various techniques, such as orthogonal or non-orthogonal transformations, covariance-based or covariance- free techniques, as will be discussed in more detail later.
  • the problem solver 120 performs sampling in a parameter space associated with the reduced base to generate one or more output model parameter sets.
  • sampling techniques may be employed to select model parameter sets in the reduced parameter space. For example, geometric or stochastic sampling may be performed. Alternatively, sampling is performed while optimizing in the reduced base. For example, Quasi-Newton methods and a moving local hyperquadric may be used.
  • geometric sampling is performed in the parameter space associated with the reduced base. Geometric sampling is especially useful where the computation of the forward model is very costly and includes sampling the coefficients of the reduced base in a bounded polytope. Different sparse sampling techniques, such as Latin hypercube sampling, Smolyak grids or random sampling techniques, may be performed in the reduced base. In addition, the sampling domain may be restricted by imposing different kinds of bounding and regularity constraints on the original model space. In another implementation, oriented sampling is performed. Oriented sampling may be performed by calculating the Jacobian on the reduced base, and using the information to orientate the sampling. Oriented sampling and geometric sampling may be applied at the same time. Oriented sampling can be viewed as a geometric sampling modality.
  • stochastic sampling is performed in the parameter space associated with the reduced base. Stochastic sampling is especially useful when the forward model is fast to compute. Stochastic sampling may be performed by using global optimization techniques, such as particle swarm optimizers, differential evolution, etc., and/or Monte Carlo methods. The search may first be performed in a prior reduced base until the valley of equivalent model parameter sets in the posterior is located. Once the valley is located, a reduced model posterior base may be used to improve the sampling efficiency.
  • the output model parameter sets are analyzed to generate one or more measures of interest that are useful for further study and assessment.
  • the output model parameter sets are analyzed to generate one or more measures of interest that are useful for further study and assessment.
  • the uncertainty assessor 122 generates uncertainty measures (e.g., indicator maps, percentile maps, E-types, etc.) of the output model parameter sets. These measures indicate, for example, the probability of occurrence of a certain feature of the model parameters in the region of low misfits. They may be used for assessing risk in various applications. For example, in geophysical applications, such measures may indicate the most likely hydrocarbon per volume, the probability of high resistance, etc. More details will be provided in the following sections.
  • uncertainty measures e.g., indicator maps, percentile maps, E-types, etc.
  • the method of reducing the base may depend on the prior knowledge (e.g., qualitative or informative) that is provided by the client (at e.g., client workstation 134).
  • the choice of model reduction approach depends on the components that are available in the input model parameter set.
  • the reduced base may be a prior base or a posterior base.
  • a prior reduced base is constructed based on an ensemble of input model parameter sets (mi, m 2 , m3, ... ), whereas a posterior reduced base is constructed based on a unique model parameter set mo having low misfit.
  • the posterior reduced base offers the advantages of higher computational efficiency and easy adaptability to the state of the art of the client's interpretation.
  • the client may not want to change the input interpretation of the problem and would like to explore uncertainties around a prescribed model parameter set.
  • the client may not have a specific interpretation in mind, and the sampling has to proceed from the beginning (i.e. prior or posterior base with optimization).
  • the reduced base may be covariance-type or co variance-free.
  • Covariance is a measure of the linear relation between random variables.
  • Examples of covariance-type reduced bases include bases that are constructed by employing a PCA- based technique (e.g., PCA or kernel PCA) using a set of prior scenarios or Singular Value Decomposition (SVD) using a unique model parameter set.
  • Govariance-free reduced bases include those constructed by using, for example, Indicator base, Discrete Cosine Transform (DCT), Discrete Wavelet Transform (DWT), Discrete Chebyshev Transform (DCHT), or Independent Component Analysis (ICA).
  • DCT Discrete Cosine Transform
  • DWT Discrete Wavelet Transform
  • DCHT Discrete Chebyshev Transform
  • ICA Independent Component Analysis
  • Covariance-free bases advantageously avoid the need for matrix diagonalization, which may be infeasible in very high dimensional spaces. It is understood that any other discrete transforms based on orthogonal polynomials (e.g., Krawtchouk, Hahn, Poisson-Charlier, etc.), may also be used to find the reduced base.
  • orthogonal polynomials e.g., Krawtchouk, Hahn, Poisson-Charlier, etc.
  • a prior reduced base may be constructed based on multiple prior scenarios. Constructing a prior reduced base is particularly useful where only qualitative information is available for generating a set of unconditional prior scenarios. This includes the situation where a set of gridded secondary attributes that contain partial information about the structure of the inverse problem solution (e.g., model parameters) is available, but the degree in which this set of attributes informs the solution is not known. The prior reduced base may then be constructed using these secondary attributes.
  • a set of gridded secondary attributes that contain partial information about the structure of the inverse problem solution (e.g., model parameters) is available, but the degree in which this set of attributes informs the solution is not known.
  • the prior reduced base may then be constructed using these secondary attributes.
  • a prior reduced base is constructed using Principal Component Analysis (PCA).
  • PCA is useful for determining orthogonal bases using different unconditional prior scenarios.
  • PCA involves a mathematical procedure that transforms a number of possibly correlated variables into a number of uncorrelated variables called principal components, which are related to the original variables by an orthogonal transformation (i.e. covariance-based).
  • PCA includes finding an orthogonal base of the experimental covariance matrix estimated with these prior scenarios, and then selecting a subset of the most important eigenvalues and associated eigenvectors that are used as a reduced model space base. The resulting transformation is such that the first principal component accounts for as much of the variability and each succeeding component accounts for as much of the remaining variability as possible. Kernel PCA methods may also be applied to calculate the experimental covariance and establish the model reduction.
  • a PCA-reduced base provides an optimal compression by gathering the commonality between prior scenarios.
  • Linear inverse theory allows computing a PCA- reduced base from the eigenvalue decomposition of matrix JF.J, JF , known as linearized Hessian.
  • the PCA-reduced base may be used to optimize and find a good model parameter set in the valley.
  • the search for the good model parameter set may be performed using a global optimization method (e.g., particle swarm optimization).
  • the geometric sampling approach may also be performed directly on the PCA-reduced base, as will be discussed later.
  • the matrix C 1 '" (JF ⁇ JF mo ) ⁇ is the pseudo-inverse (see symbol (») ⁇ ) of the linearized Hessian JF' JF , and it is related to the posterior covariance matrix in the model parameter set m 0 .
  • Matrix C 1 '" is a symmetric and definite positive matrix, and thus admits orthogonal diagonalization as follows:
  • Target oriented bases are built from prior scenarios and learning which terms are affected the most by the nonlinearities of the forward operator. This generally includes decomposing (or projecting) the different prior scenarios into the prior reduced base and learning how to modify these coefficients in order to minimize the data misfit, taking into account the misfits of the prior scenarios that have been previously calculated.
  • the learning algorithms may be important when the forward model is very costly, since these prior scenarios have been already pre- computed to determine their misfit.
  • stochastic optimization techniques may be used instead.
  • Geometric sampling techniques may also be used for costly forward solves. Nevertheless, the sampling process may be inefficient if the base is constructed using only prior knowledge.
  • a set of prior scenarios are first generated using, for example, conditional simulation techniques.
  • the misfit values for the prior scenarios are then computed.
  • the prior scenarios are projected onto the reduced base to determine their coefficients.
  • learning methods are used to unravel the structure of these coefficients in order to fit the observed data. Sampling may then be performed in the region of equivalence using the target oriented bases and sampling techniques that will be described later.
  • the reduced base there are several ways of determining the reduced base, depending on the type of prior information or knowledge that is at our disposal or the cost of individual forward solves.
  • a set of unconditional prior scenarios may be generated by such qualitative information.
  • the prior reduced base may be built using well log information and secondary geophysical attributes.
  • a PCA base provides a reduced base that gathers the commonality between these priors, thereby providing an optimal model reduction.
  • the present framework accommodates the use of secondary attributes where data of the primary attributes of interest is not available or missing.
  • the prior information may be, for example, a set of gridded secondary attributes which are correlated to the model parameter set of interest, and therefore provides partial information about the structure of the solution to the predetermined problem. For example, assuming that the interest lies in determining the uncertainty of meteorological data, but precipitation data in a certain region (i.e. primary attribute or model parameter set) may not be available. Nevertheless, the PCA base (or target oriented base) may still be constructed using a set of different secondary attributes, such as the altitude of the region and other suitable digitized attributes (e.g., insolation, slope, etc.).
  • a more informative secondary attribute may be available.
  • the informative secondary attribute may be completely sampled in a grid. For instance, the porosity in a reservoir is well correlated to the seismic acoustic impedance and may be considered as an informative secondary attribute.
  • Such secondary attribute provides information about the approximate structure of the solution of the predetermined problem (or model parameter sets of interest).
  • the hard model information may be precipitation levels at certain meteorological locations, wherein the interest lies in determining the reduced base of the precipitation (or model parameter set) of a particular other region.
  • the set of prior scenarios may be built using image-guided interpolation and/or condition simulation techniques based on the hard model information and/or secondary attribute (e.g., altitude).
  • p terms of the s k base may be selected in order to fit in, for example, the least square sense, the hard model information m ( ⁇ a ), - l, N p .
  • This spatial trend may be used to calculate the misfit in the hard model positions.
  • a set of conditional simulations of this residual field may be generated.
  • PCA-based technique may be employed to expand the high frequency content of the model.
  • telescopic sampling may be performed using t as a trend model and
  • a posterior reduced base is constructed based on a unique posterior model parameter set mo.
  • the posterior reduced base may be orthogonal-based (e.g., SVD, DCT, DWT, DCHT, ICA, etc.) or non-orthogonal based (e.g., indicator bases, etc.).
  • FIG. 3a shows an exemplary method 204 of constructing an orthogonal reduced base.
  • an orthogonal reduced base is constructed by applying one or more orthogonal transformations to the posterior model parameter set mo.
  • a unique posterior model parameter set mo is received.
  • the model parameter set mo may be generated as previously described.
  • the model parameter set mo may be a 2-D image 305 that serves as a candidate solution to an inverse problem. It is understood that although the present framework is described with reference to 2-D images, it is also applicable to parameter sets associated with one, three or other dimensionalities where separability exists or where the reduction method allows for it.
  • orthogonal transformation matrices are defined. In one
  • U and V are orthogonal transformation matrices related by ⁇ , that is diagonal by blocks with non- negative real numbers.
  • DCT digital coherence tomography
  • DWT dense coherence to DCHT
  • ⁇ is not diagonal by blocks. Instead, its non-null coefficients are concentrated on its first rows and columns.
  • the orthogonal transformation matrices (or bases) U and V may be determined from the matrices m 0 m 0 7 , m 0 respectively.
  • orthogonal transformation matrices are possible, depending on the model reduction method used.
  • model reduction methods include, for example, Singular Value Decomposition (SVD), Discrete Cosine Transform (DCT), Discrete Wavelet Transform (DWT) using different families of wavelets, Discrete Chebyshev Transform (DCHT), as will biTdescribed later.
  • any other discrete transforms based on orthogonal polynomials e.g., Krawtchouk, Hahn, Poisson-Charlier, etc.
  • non-orthogonal model base reduction may also be performed using indicator or level set bases. More details will be provided in the following sections.
  • the model parameter set m 0 is compressed by projecting onto one or both orthogonal matrices (U and/or V) to generate compressed terms (Ui mage and/or Vimage) as follows:
  • the model parameter set m 0 may be compressed horizontally or vertically respectively, as illustrated by exemplary compressed images 31 1a and 311b.
  • the compressed term U ma ge is the result of uncorrelating and compressing the first rows r d of U, corresponding to most of its variability.
  • the compressed term is the result of uncorrelating and compressing the first columns c d of .
  • V image or ⁇ image depends on, for example, the nature of the model parameter set m 0 .
  • U may provide a more effective compression than V, where the model parameter set m 0 is horizontally stratified.
  • ⁇ image mav provide a more effective compression than U image where the model parameter set m 0 is vertically stratified.
  • a thresholding matrix ⁇ may be defined as a zero matrix with the dimensions of U, but containing only its z ' -th row:
  • the thresholding matrix T ( may be defined as a zero matrix with the dimensions of V, but containing only its z ' -th column:
  • the reduced base may be derived using the thresholding matrix ⁇ , ⁇ by projecting the compressed terms back onto the original canonical base associated with the posterior model parameter set m 0 .
  • This exemplary model reduction using the orthogonal matrix V is particularly useful when U is hard to compute (e.g., portrait images in overdetermined systems).
  • C j has the same role as r d in the case of V image .
  • the same transformation may also be performed row-wise as follows:
  • the reduced base generated in accordance with the present framework offers many useful properties that can be exploited.
  • One such advantageous property for example, is that exact model compression may be performed.
  • Exact model compression involves transforming back all the higher frequencies vectors of the U, (or V ⁇ mfl ) into a single base vector that gathers all the high frequency terms of the orthogonal decomposition. For instance, in the case of horizontal compression using U moj?e , the transformation is as follows:
  • Group compression includes compressing groups of rows or columns of the uncorrected image represented by the threshold matrices.
  • the transformation for group compression may be as follows:
  • Group compression is particularly useful when discrete wavelet transform (DWT) is applied, because the base terms have a compact support that varies with the wavelet family that has been used. It is also useful when Chebyshev transform is applied, because Chebyshev polynomials of degree higher than fifty are unstable. In such case, the base is calculated by transforming the image by blocks.
  • DWT discrete wavelet transform
  • Chebyshev transform is applied, because Chebyshev polynomials of degree higher than fifty are unstable. In such case, the base is calculated by transforming the image by blocks.
  • FIG. 3b shows images 330a-d representing exemplary reduced base terms with increasing spatial frequencies. As the spatial frequency of the higher terms of the reduced base increases, the reduced base terms 330a- d expand details of the image corresponding to more uncertain parts of the model space. Although a certain number of base terms may be used to project the image (or model parameter set), all of them may be at disposal and may be used to analyze the uncertainties related to the presence of these high frequency terms. This advantageous property is useful in the telescopic geometric sampling approach, as will be described in more detail later.
  • Singular Value Decomposition is employed to define the orthogonal transformation matrices in step 306 of FIG. 3a.
  • SVD represents the factorization of any rectangular matrix G e (s, w) in the following form:
  • U and V are orthogonal matrices that provide orthonormal bases for sets of real numbers (R s and R") in s- and ⁇ -dimensional spaces respectively.
  • ⁇ is a s-by-n box- diagonal matrix with non-negative real numbers on the diagonal, called the SVD of the matrix (i.e. image) G .
  • the base U comes from the eigenvalue decomposition (e.g., PCA) of the row correlation matrix GG T
  • the base V may be calculated from the eigenvalue decomposition (e.g., PCA) of the column correlation matrix G T G .
  • the row and column correlation matrices GG T and G T G may be calculated as follows:
  • the model parameter set m 0 is a 2-D image, such as the conductivity or the P-velocity field in a 2-D electromagnetic or seismic tomography problem.
  • FIG. 4 illustrates an exemplary application of SVD to a 2-D representation of a geological layer (or model) to achieve model reduction.
  • the original image 402 representing m 0 may be reconstructed into image 404.
  • the reconstruction is based on the first 5 principal modes of the image 402 determined using SVD.
  • Images 406a-d show the first 4 terms of the reduced base.
  • the image m 0 may be projected onto any or both bases, as described previously with respect to step 310 shown in FIG. 3a:
  • a matrix, T may be defined as a zero matrix with the dimensions of U, but containing only its z-th row:
  • the reduced base may be derived by projecting the first r d threshold matrices back onto the original canonical base as follows:
  • SVD achieves optimum compression, but it may not be feasible to apply on a whole image with large number of pixels. Some free-diagonalization methods may be required to achieve the compression and find the reduced base.
  • SVD may be applied to both 2-D and 3-D images. In one implementation, higher order SVD may be used to expand this methodology to 3-D.
  • the 3-D image may be transformed into a new 2-D image by reshaping its different layers. The reduced base may be determined based on this new 2-D image and transformed back to the original 3-D space.
  • the U base may be generated by diagonalizing m 0 m 0 T following the smallest dimension.
  • V T image m " ⁇ V T ⁇
  • DCT discrete cosine transform
  • DFT discrete Fourier transform
  • DCT is a separable transformation that may also be defined in higher dimensions.
  • FIG. 5 shows the exemplary application of DCT on an image 502.
  • Image 504 may be reconstructed by using the first 5 principal modes of the image determined by using DCT.
  • Images 506 a-d show the 4 first terms of the reduced base.
  • DCT is written in matrix form as follows:
  • D is the two-dimensional cosine transform of the m 0 image
  • L and R are two orthogonal matrices, which may be defined as follows:
  • discrete wavelet transform is employed to define the orthogonal transformation matrices in step 306 of FIG. 3a.
  • FIG. 6 shows an exemplary application of DWT on an original image 602, which represents a geological layer.
  • the exemplary image 604 may be reconstructed by using the first 4 principal modes (606a-d) of the image 602, which are determined by DWT using the Daubechies family.
  • wavelet families e.g., Symlets, Coiflets, etc.
  • the original image (or model parameter set) is represented in the form m 0 e M s,n .
  • the image mo may be decomposed into orthogonal wavelet transformation matrices V w and V render, as follows: (47)
  • the image, m 0 may be projected onto one or the other (Untended,and V w ) for horizontal or vertical compression respectively, as follows:
  • G G where H represents a low pass or averaging portion and G represents the high pass or differencing portion of the matrices Wi and WR.
  • H represents a low pass or averaging portion
  • G represents the high pass or differencing portion of the matrices Wi and WR.
  • V is the vertical difference
  • H is the horizontal difference
  • D is the diagonal difference
  • U image may be expressed as follows:
  • V mi3 ⁇ 4e
  • V m *" ⁇ W " « T (Hm 0 T ) T
  • (Gm 0 T ) T ⁇ [ra 0 H T
  • Wavelets have compact support to allow localization of singularities in the space (time)-frequency domain. Due to the compact support of the wavelets, the wavelet reduced bases may also have a compact support. To determine the reduced base, certain groups of spatial frequencies may be transformed together back to the original space.
  • DCHT Discrete Chebyshev Transform
  • the orthogonal transformation matrices in step 306 of FIG. 3a are determined using discrete Chebyshev transform (DCHT).
  • DCHT discrete Chebyshev transform
  • FIG. 7 shows an exemplary application of DCHT to the image 702, which serves as a 2-D input model parameter set mo.
  • the exemplary image 704 may be reconstructed by applying DCHT to determine the first 5 principal modes of the image 702.
  • Image 706a-d show the first 4 terms of the DCHT reduced base.
  • DCHT is a separable transformation that may also be defined in one dimension (1-D), two-dimensions (2-D) or higher dimensionalities.
  • t u (/ ' ) , t v [j) are 1-D discrete scaled Chebyshev polynomials as defined below.
  • N usually the matrix size
  • x in the range [0, N - 1]
  • the model reduction may be performed by constructing two orthogonal matrices, U s and V n , such that:
  • the reduced base may be calculated by transforming the image by blocks.
  • ICA independent component analysis
  • FIG. 8 shows the exemplary application of ICA on the original image 801. Reconstruction may be performed by applying ICA on the first 4 independent principal base terms (802a-d).
  • ICA Independent component analysis
  • both the linear transform A and the independent components Ic are determined by maximizing the non-Gaussian character of the independent components.
  • FIG. 9 illustrates an exemplary alternative method 204 for computing an orthogonal reduced base by using spatial indicators.
  • the computation of this base is not based on orthogonal transformations, as previously described with respect to FIG. 3a.
  • the method 204 includes generating the indicator base (or level set) of an image Image e M(s, n) .
  • a posterior model parameter set is received.
  • the posterior model parameter set is generated using the previously described methods.
  • the model parameter set may be a 2-D image ("Image"), such as geological layer image 903. Other dimensionalities (e.g., 3-D) are also useful.
  • the percentile distribution of the model parameter set (or Image) is computed.
  • the percentile distribution describes the statistical distribution of a random variable (e.g., pixel values) associated with the image.
  • the model percentile distribution may be represented by, for example, an empirical cumulative distribution function (CDF), such as graph 905 shown in FIG. 9.
  • CDF empirical cumulative distribution function
  • the mean ⁇ may be computed globally, column-wise and/or row-wise, depending on the desired amount of variability in the image that needs to be maintained. In such case, the percentile distribution may be based on the residual image. For each pixel of the residual image R, an indicator cut-off variable is defined for a certain cut-off threshold t:
  • the spatial indicator base is determined using the spatial indicator terms ⁇ , , ⁇ . ,..., ⁇ ⁇ .
  • two consecutive terms of this set are subtracted to produce the reduced orthogonal base
  • the residual image R obtained after subtracting the mean is projected onto this base by, for example, the least-squares method:
  • M is the matrix lb, ,Ib 2 ,...,Ib ? J built with the vectors of the indicator base.
  • the last term or residual Ires may be computed so as to make the reconstruction exact:
  • Ires Imager - ⁇ - ⁇ ⁇ , ⁇ , (63) where q is the number of percentiles that have been selected and gives the dimension of the reduced space into which the model parameter set is projected. If the residual Ires is taken into account, the dimension of the reduced space will be q+1.
  • FIG. 10 illustrates the exemplary application of the spatial indicator-based method 204 on the original image 903.
  • Image 1004 may be reconstructed by applying the spatial indicator-based method 204 on the first 5 principal modes of the original image 903.
  • Images 1006a-d show the first 4 terms of the spatial indicator reduced base.
  • the reduced base facilitates the computation of the Jacobian and Hessian matrices, which were previously almost impossible or infeasible to compute in high dimensional spaces.
  • the Jacobian matrix represents the first-order derivatives of the cost function in the reduced base, while the Hessian matrix is a square metric of the second-order partial derivatives. Both matrices are useful in the sampling and/or optimization steps, which will be described later.
  • the Jacobian for example, may be used to provide the local direction of the valley and orientate the sampling.
  • the Hessian provides curvature information of the valley, and the kind of hyperquadric that may be locally fitted to the cost function.
  • the Hessian also provides information about the presence of saddle points that serve to connect different basins of potential solutions.
  • the q coefficients may be perturbed around the point (l,...,l)j v v j as follows: d F F((l,-,l, +g ,-,l ? ) B J-F((l span-,l ,-,l,) g
  • the Jacobian matrix JF ⁇ ⁇ U ⁇ s s) ⁇ (s g) ( ⁇ >9) ) may be generated by computing q additional forward solves for all the data points.
  • the Jacobian matrix may have a size (s,q), and if s>q and its rank is s, its singular value decomposition may have non-null singular values that will give us the extent of the equivalences on the reduced space.
  • the q singular vectors contained in V ( r q q) provide the local orientation of the valley in the reduced space.
  • Equation (64) means that the problem becomes less ill-posed in the reduced space. This is because s will become greater than q eventually, and all the singular values of the Jacobian will be non-zero, thereby creating an over-determined system. This provides a criterion to determine, in practice, the maximum number of terms in the base that are going to be informed by the data (s).
  • Relationship (68) may also be used to express the Jacobian in the reduced base once it is known in the pixel base. This may be important because the Jacobian in the pixel base may also be calculated using adjoint methods. Furthermore, the Hessians of the component functions f k in both bases are orthogonally similar, as shown by the following relationship:
  • sampling may be performed in the reduced- dimensional model parameter space associated with the reduced base using various sampling methods. These sampling methods include geometric sampling, stochastic sampling, oriented sampling, and sampling while optimizing using the local hyperquadric, as will be described in the following sections. It is understood that other types of sampling methods may also be applied.
  • Geometric hyperprism sampling may be performed to generate an output model parameter set.
  • it is feasible to apply geometric sampling even if the original model parameter space has many dimensions and the forward model is very costly to compute, because the sampling and the forward evaluation steps are advantageously decoupled.
  • Sampling may include different sparse sampling techniques, such as Latin hypercube sampling, Smolyak grids, random sampling or other schemes.
  • the sampling domain is restricted by imposing different kinds of bounding and regularity constraints on the original model parameter space.
  • oriented sampling may be performed by computing the Jacobian on the reduced base and using this information to orientate the sampling.
  • geometric sampling includes bounded constrained sampling.
  • this approach involves imposing lower and upper bounds (1, u) in the original model parameter space.
  • a search hyperprism may be circumscribed in the reduced base.
  • the vertices of the search hyperprism represent the extrema of a region of allowable coefficients a k in the reduced base.
  • the search hyperprism may be circumscribed by directly solving a vertex enumeration problem, or by other computationally feasible methods, such as axes projection, optimization or sparse random sampling.
  • FIG. 1 1 shows an exemplary bounded constrained sampling method 206.
  • the reduced base is received.
  • the reduced base may be a posterior or a prior reduced base.
  • the reduced base is associated with a set of orthogonal vectors Vk.
  • the reduced base may be generated by, for example, step 204 previously described with respect to FIGS. 2-10.
  • the model parameter set m is bounded by lower and upper bounds (I, u) in the original model parameter space.
  • the lower and upper bounds (1, u) may be defined based on the prior knowledge of the smallest and largest allowable or expected properties.
  • ⁇ is the model mean
  • a k are real valued coefficients
  • v k are the vectors of the reduced base.
  • the bounding constraints associated with the coefficients a k may be implemented by circumscribing a search hyperprism in the reduced base.
  • the search hyperprism may be circumscribed by directly solving a vertex enumeration problem, which may include solving the linear upper and lower constrained problem expressed by the above-mentioned equation (71).
  • the solution of the vertex enumeration problem generates a ⁇ -dimensional convex polytope in R ? , which defines the sampling search space.
  • the circumscribed hyperprism may also be determined without explicitly solving the vertex enumeration problem.
  • the hyperprism is circumscribed by axes projection or coordinate search on each individual coefficient a k .
  • Axes projection may be performed by projecting upper and lower bounds (1, u) into the reduced space, as follows: mJ L ⁇ L) ⁇ a , (72)
  • sampling is performed within the bounded parameter space (e.g., circumscribed hyperprism and/or polytope) to generate the coefficients a k .
  • Various sampling procedures such as the Latin hypercube sampling, sparse collocation grid sampling (e.g., Smolyak grids with Clenshaw-Curtis or Gaussian abscissas) or even random samplers, may be used for this purpose.
  • Very sparse sampling methods are particularly useful in computationally expensive forward problems. More exhaustive sampling procedures (e.g., uniform sampling) may also be used where there are few dimensions and fast forward solvers.
  • Sampling may also be performed without first circumscribing the hyperprism and/or polytope.
  • the coefficients a k are determined by sparse random sampling while imposing, at the same time, the bounding constraints (71). Sampling either with or without first circumscribing the hyperprism and/or polytope effectively decouples the sampling procedure from the forward evaluations.
  • the resulting samples i.e. coefficients a k
  • classification may be performed using clustering techniques in the reduced space before any forward evaluation. This procedure is particularly useful for very computationally expensive forward or inverse problems.
  • nti denotes the model parameter sets and c ; denotes the model thresholds.
  • constraints may also be added, such as in the oriented sampling approach, in which sampling is constrained to increase its efficiency.
  • the results of the geometric sampling are returned.
  • the results may include the generated samples (i.e. coefficients ⁇ 3 ⁇ 4) and the associated classification and probabilities in the original base.
  • the sampling results after likelihood evaluation may be further processed by step 208 in FIG. 2, as will be described later.
  • geometric sampling approaches may also be used. These geometric sampling variants may be defined by how the reduced base vectors v k and the mean term ⁇ are chosen. For example, minimum norm geometric sampling may be applied. In such case, the input model parameter set mo is computed using the Moore- Penrose pseudo inverse of the linearized problem. Accordingly, this model parameter set does not have components in the local null-space of the forward operator F. The uncertainties are determined by sampling in these directions.
  • the minimum norm geometric sampling method includes determining the minimum norm solution of the linearized problem, m ⁇ , and applying the model reduction techniques to and determining the terms of the reduced base.
  • FIG. 12 shows an exemplary telescopic sampling method 206, which is another implementation of a geometric sampling-based method. Telescopic sampling is possible in the present framework, because of the orthogonal character of the reduced base that allows high frequencies may be added in the sampling as needed.
  • the prior reduced base is received.
  • the reduced base is associated with L + Q base terms.
  • the reduced base terms may be computed by, for example, step
  • the mean model term ⁇ is determined using the first L base terms.
  • the mean model ⁇ ⁇ is the expansion of a "good" model parameter set using a certain
  • the present framework advantageously enables the use of other model parameterizations (e.g., logarithmic sampling) that facilitate better conditioning of the sampling procedure.
  • model parameterizations e.g., logarithmic sampling
  • the logarithmic parameterization helps to linearize the search.
  • trade-offs between model parameter sets in a lot of cases are products and quotients, and these relationships become linear in the logarithmic space.
  • the logarithmic reduced geometric sampling includes applying the model reduction technique to the logarithm of the model and performing the search on the logarithmic space.
  • the bounds are in the logarithm space, and the Jacobian and the Hessian in the reduced space accounts for the logarithmic transformation using the chain rule, as stated above.
  • Oriented sampling may be performed alternatively or in combination with geometric sampling (at 206), and may also be viewed as a geometric sampling modality.
  • oriented sampling is performed by computing the Jacobian of the forward operator F in the reduced base and using this information to orientate the sampling.
  • FIG. 13 shows a 2-dimensional plot 1302 of the valleys of equivalent model parameter sets.
  • the Jacobian of the forward operator F in the reduced base provides the local directions 1304 of the valley 1306.
  • oriented sampling involves solving the quadratic problem of finding the coefficients a k expressed in the following exemplary relationship:
  • Relationship (71 ) may be added as a constraint in the search if needed.
  • the following approximation may be directly applied:
  • oriented geometric sampling is performed by adding the Jacobian-based constraint (78 or 81) to the geometric sampling approach expressed by, for example, relationship (71). The Jacobian-based constraint serves to restrict the polytope in the reduced space.
  • oriented geometric sampling is performed by directly sampling on the J k base imposing I ⁇ a k J k ⁇ u .
  • the logarithmic parameterization facilitates the linearization of the region of nonlinear equivalence, thereby enhancing the efficiency of the oriented sampling procedure.
  • J k terms provide the local orientation of the valley
  • global bounds may be imposed to facilitate a global search along the valley.
  • Stochastic sampling may also be performed in the reduced base to generate an output model parameter sets. Stochastic sampling is particularly useful when the forward model is fast to compute. It may involve the use of optimization techniques such as Particle Swarm optimization (PSO), differential evolution, Monte Carlo methods, or a combination thereof.
  • PSO Particle Swarm optimization
  • Monte Carlo methods or a combination thereof.
  • stochastic sampling is performed in a posterior reduced base.
  • FIG. 14a shows an exemplary method 204 of computing the posterior reduced base
  • FIG. 14b shows an exemplary method 206 of performing stochastic sampling in the posterior reduced base.
  • prior information is received.
  • Such prior information may include, for example, prior unconditioned scenarios (or different model parameters sets) that do not fit the observed data.
  • the prior information is used to build a prior reduced base.
  • the prior reduced base comprises a PCA base.
  • the prior reduced base may be computed based on an ensemble of prior model parameter sets.
  • optimization is performed in the prior reduced base until a valley of equivalent model parameter sets (or low misfit region) in the posterior is located.
  • the optimization may be performed with any kind of regularization.
  • particle swarm optimization is performed.
  • PSO may include using members of the PSO family with the highest convergence rate (e.g., regressive- regressive-PSO) until a good model parameter set m b in the valley of equivalent model sets is found.
  • Other PSO family members or other optimization methods e.g., differential evolution may also be used.
  • the posterior reduced base is computed.
  • the posterior reduced base is determined by applying model reduction techniques using, as
  • the model reduction techniques include, for example, SVD,
  • FIG. 14b shows an exemplary stochastic method 206 that may be performed in the posterior reduced base.
  • the posterior reduced base is received. By switching to the posterior reduced base, more efficient sampling may be performed.
  • sampling is performed in the posterior reduced base using a global optimization technique.
  • exemplary global optimization techniques include, for example, PSO (e.g., explorative versions such as centered-progressive-PSO), differential evolution (or sampling while optimization), or Monte Carlo methods (or importance sampling). Other types of global optimization methods may also be used.
  • PSO e.g., explorative versions such as centered-progressive-PSO
  • differential evolution or sampling while optimization
  • Monte Carlo methods or importance sampling.
  • Other types of global optimization methods may also be used.
  • the sampling and the forward evaluations are coupled.
  • Some PSO family members e.g., regressive- regressive-PSO, centered-centered-PSO and generalized PSO
  • PSO variants such as centered-progressive (CP)- PSO, progressive-progressive (PP)-PSO and regressive-regressive (RR)-PSO, have also achieved good exploration of the search space.
  • posterior analysis of the likelihood of collected samples is performed. For example, posterior statistics such as marginals, e-types, indicator maps, etc, may be generated based on the ensemble of output equivalent model parameter sets.
  • the results of the stochastic sampling are returned for further processing by, for example, step 208 of FIG. 2.
  • Sampling may also be performed while optimizing in the reduced base. More particularly, local optimization and hyperquadric sampling in the reduced base may be iteratively performed until a stop criterion is fulfilled. Local optimization methods, such as Quasi-Newton methods (e.g., Gauss-Newton methods), may be performed to compute a moving local hyperquadric that circumscribes the sampling search space, by finding the minimum of the function
  • Quasi-Newton methods e.g., Gauss-Newton methods
  • the Gauss-Newton method is used to perform the local optimization.
  • the Gauss-Newton method may be used to minimize a proxy of the function
  • J3 ⁇ 4 i JF. i (m - m 0 ) JF: i (F(m 0 ) -d).
  • JF' JF m may be referred to as the linearized Hessian, because the second order terms have been removed in the first order approximation of F (m) .
  • FIG. 15 illustrates an exemplary method 206 of sampling while optimizing in the reduced base in further detail.
  • a reduced base is received.
  • the reduced base may be a posterior or a prior reduced base.
  • a posterior reduced base for example, may be generated based on a given input model parameter set m 0 , as follows:
  • the Jacobian of the cost function C(m) is computed in the reduced base.
  • the Hessian HC ⁇ a j may be defined as follows:
  • a local hyperquadric is computed (or circumscribed) by using the Hessian and sampling is performed within the hyperquadric by circumscribing an oriented hyperprism.
  • the axes of the hyperquadric, rotated to the V base, are the inverse of the eigenvalues of the Hessian.
  • the hyperquadric may be regularized by adding to the smallest eigenvalues a tolerance that is calculated in accordance with the ratio to the biggest eigenvalue.
  • the hyperquadric may be sampled using, for example, sparse grids or any other sampling method. Calculations are performed in the reduced base.
  • steps 1504- 1508 are repeated. In particular, the center of the next hyperquadric is processed. Steps 1504-1508 are iteratively repeated until a predetermined stop criterion is fulfilled.
  • Exemplary stop criteria are associated with, for example, error tolerance, the number of function evaluations, the number of sampled points, and so forth.
  • the method 206 checks the value of the cost function in the centre of the new
  • the reduced base and Jacobian of the cost function at the center of the hyperquadric are recalculated.
  • the Hessian may also be recalculated, in the reduced base, at the new center.
  • the Hessian may be approximated using a numerical approximation method, such as a Quasi-Newton method (e.g.,Davidon- Fletcher-Powell method (DFP), Broyden-Fletcher-Goldfarb-Shanno (BFGS) method, limited memory BFGS method (L-BFGS) or the Symmetric Rank 1 method (SRI)). All these methods have scalable formulae to calculate the inverse of the Hessian, which may be used to compute the center of the local hyperquadric.
  • DFP Davidon- Fletcher-Powell method
  • BFGS Broyden-Fletcher-Goldfarb-Shanno
  • L-BFGS limited memory BFGS method
  • SRI Symmetric Rank 1 method
  • the results of the sampling are output for further processing by, for example, step 208 of FIG. 2.
  • FIG. 16 shows a 2-D plot 1602 illustrating this method for the Rosenbrock function, which is a synthetic benchmark function that has a valley shape. It is also known as the "banana” function. For sake of simplicity, the hyperprisms are not drawn.
  • the output model parameter sets are analyzed at 208 to generate one or more measures of interest.
  • the output model parameter sets are analyzed at 208 to generate one or more measures of interest.
  • one or more uncertainty measures or probability maps are computed.
  • exemplary uncertainty measures include indicator maps, percentile maps or E-type values. It should be understood that other measures of interest may also be identified for further study.
  • Indicator maps may be generated by computing, at each pixel of the model parameter set (m), the probability that the model parameter set is less than a certain cut-off c, as follows:
  • percentile maps are generated by computing, at each pixel of the model parameter set (m), the / ⁇ -percentile of the distribution of the values at this pixel deduced from the ensemble empirical distribution, as follows:
  • the 50% percentile is the median, the 25% percentile is the lower-quartile, and the 75% percentile is the upper-quartile.
  • Lower and upper spatial percentiles may be used to detect percentile anomalies.
  • E-type values are generated. E-types are the mean values (or the median) in each pixel deduced from the ensemble. Standard deviations and inter-quartile ranges may also be determined and interpreted as a measure of model variability (or uncertainty).
  • Exploration geophysics involves a set of techniques that map the subsurface structure of a region to elucidate the underlying structures, spatial distribution of rock units, and to assess the likelihood of ore deposits or hydrocarbon accumulations. These techniques include, but are not limited to, seismic reflection techniques, gravimetric, magnetic or electromagnetic methods, electrical methods and tomography techniques. [00188] More particularly, seismic reflection techniques may be used to map the subsurface distribution and structure of terrains so as to delineate potential hydrocarbon accumulations. In a process known as seismic inversion, the seismic reflection data can be transformed into a quantitative rock property description of a reservoir.
  • Gravimetric, magnetic and electromagnetic methods are also very important in mineral exploration and in facilitating the delineation of the hydrocarbon reservoirs in oil and gas exploration.
  • magneto-telluric and controlled-source electromagnetic (CSEM) data can be used to detect resistivity changes, complementing seismic data to detect hydrocarbons.
  • electrical methods such as electrical resistivity tomography, induced polarization and spontaneous potential methods can be employed in geothermal studies, hydrogeology and environment.
  • Reservoir simulation is an area of reservoir engineering in which computer models are used to predict the flow of fluids (e.g., oil, water, and gas) through porous media.
  • Reservoir simulation models are used by oil and gas companies in the development of new fields, and to make production forecasts for investment decisionmaking for developed fields. Such models may be used to identify the number of wells required for exploiting a new field and to manage the way that the reservoir is going to be produced to improve oil recovery: hydraulic fracturing, deviated wells, water and gas reinjection, etc.
  • Other reservoir simulation applications include solving history matching inverse problems, which try to improve the knowledge of the reservoir properties by reproducing its past behavior (e.g., historical production or pressures).
  • model parameter set (or a family of model parameter sets) has been historically matched, it can be used to simulate future reservoir behaviour with a higher degree of confidence.
  • Other applications include well placement and net present value optimizations.
  • the model parameters include, for example, permeability, porosity or saturation. Similar applications may also be found in the hydrogeology field.
  • the present framework may be applied in the fields of biomedical applications, financial analysis, and meteorology predictions.
  • the present framework may be used to determine the uncertainty associated with data acquired by various imaging modalities, such as computed tomography (CT), electrical impedance tomography, optical tomography, magnetic resonance imaging, etc.
  • CT computed tomography
  • the present framework may be applied to estimate the uncertainly of a diagnostic conclusion determined by a predictive model.
  • Other biomedical applications include probabilistic image segmentation in medical imaging and diagnosis, electromagnetic inverse scattering problem for cancer detection, classification and prediction of biomedical microarrays (e.g., proteins, cellular structures), data mining, biomolecular structural predictions, drug design, and so forth.
  • the present framework may be used to predict financial markets, optimize portfolios, assess and manage risk, and generally to analyze investments by simulating the various sources of uncertainty affecting their value, and then determining their average value over the range of resultant outcomes.
  • the present framework may be used in weather forecasting and numerical weather prediction, contamination dispersion and air quality models, etc.

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Probability & Statistics with Applications (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Algebra (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Aiming, Guidance, Guns With A Light Source, Armor, Camouflage, And Targets (AREA)
  • Compression, Expansion, Code Conversion, And Decoders (AREA)

Abstract

L'invention porte sur un cadre d'analyse de données dans un espace à grand nombre de dimensions. Selon un mode de réalisation, des données observées et au moins un ensemble de paramètres de modèle d'entrée sont reçus. L'ensemble de paramètres de modèle d'entrée sert de solution candidate d'un problème prédéfini (par exemple, un problème d'inversion ou d'optimisation) et est lié aux données observées par l'intermédiaire d'un modèle. Afin d'obtenir une efficacité de calcul améliorée, une base réduite ayant un plus petit nombre de dimensions est déterminée sur la base de l'ensemble de paramètres de modèle d'entrée. La base réduite est associée à un ensemble de coefficients, qui représente les coordonnées de n'importe quel ensemble de paramètres de modèle dans la base réduite. Un échantillonnage est réalisé dans la base réduite afin de générer un ensemble de paramètres de modèle de sortie dans la base réduite. L'ensemble de paramètres de modèle de sortie est compatible avec l'ensemble de paramètres de modèle d'entrée et est en adéquation avec les données observées, par l'intermédiaire du modèle, dans les limites d'un seuil prédéterminé.
PCT/SG2011/000207 2010-06-14 2011-06-09 Analyse de données à grand nombre de dimensions WO2011159255A2 (fr)

Applications Claiming Priority (6)

Application Number Priority Date Filing Date Title
US35468510P 2010-06-14 2010-06-14
US61/354,685 2010-06-14
US13/041,423 2011-03-06
US13/041,423 US8688616B2 (en) 2010-06-14 2011-03-06 High-dimensional data analysis
US201161494403P 2011-06-07 2011-06-07
US61/494,403 2011-06-07

Publications (2)

Publication Number Publication Date
WO2011159255A2 true WO2011159255A2 (fr) 2011-12-22
WO2011159255A3 WO2011159255A3 (fr) 2012-11-29

Family

ID=45348778

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/SG2011/000207 WO2011159255A2 (fr) 2010-06-14 2011-06-09 Analyse de données à grand nombre de dimensions

Country Status (1)

Country Link
WO (1) WO2011159255A2 (fr)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9746577B2 (en) 2013-06-03 2017-08-29 Exxonmobil Upstream Research Company Uncertainty estimation of subsurface resistivity solutions
CN107728629A (zh) * 2017-09-19 2018-02-23 富平县韦加无人机科技有限公司 无人机磁异常检测系统及方法
CN109997131A (zh) * 2016-10-26 2019-07-09 谷歌有限责任公司 用于基于内核的机器学习的结构化正交随机特征
CN110990757A (zh) * 2019-12-05 2020-04-10 杭州电子科技大学 利用无相位数据解决高度非线性电磁逆散射问题的方法
CN112559308A (zh) * 2020-12-11 2021-03-26 广东电力通信科技有限公司 一种基于统计模型的根告警分析方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4142311A (en) 1976-09-13 1979-03-06 Lane Ronald S Permanent calendar

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011115921A2 (fr) * 2010-03-19 2011-09-22 Schlumberger Canada Limited Estimation d'incertitude pour problèmes inverses non linéaires à grande échelle utilisant échantillonnage géométrique et compression de modèle sans covariance

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4142311A (en) 1976-09-13 1979-03-06 Lane Ronald S Permanent calendar

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9746577B2 (en) 2013-06-03 2017-08-29 Exxonmobil Upstream Research Company Uncertainty estimation of subsurface resistivity solutions
CN109997131A (zh) * 2016-10-26 2019-07-09 谷歌有限责任公司 用于基于内核的机器学习的结构化正交随机特征
CN109997131B (zh) * 2016-10-26 2024-01-12 谷歌有限责任公司 为基于内核的机器学习系统生成输入的方法和电子装置
CN107728629A (zh) * 2017-09-19 2018-02-23 富平县韦加无人机科技有限公司 无人机磁异常检测系统及方法
CN110990757A (zh) * 2019-12-05 2020-04-10 杭州电子科技大学 利用无相位数据解决高度非线性电磁逆散射问题的方法
CN110990757B (zh) * 2019-12-05 2023-09-19 杭州电子科技大学 利用无相位数据解决高度非线性电磁逆散射问题的方法
CN112559308A (zh) * 2020-12-11 2021-03-26 广东电力通信科技有限公司 一种基于统计模型的根告警分析方法
CN112559308B (zh) * 2020-12-11 2023-02-28 广东电力通信科技有限公司 一种基于统计模型的根告警分析方法

Also Published As

Publication number Publication date
WO2011159255A3 (fr) 2012-11-29

Similar Documents

Publication Publication Date Title
US8688616B2 (en) High-dimensional data analysis
US9619590B2 (en) Uncertainty estimation for large-scale nonlinear inverse problems using geometric sampling and covariance-free model compression
Fernandez Martinez et al. On the topography of the cost functional in linear and nonlinear inverse problems
Wang et al. A segmentation approach for stochastic geological modeling using hidden Markov random fields
Kyriakidis et al. A geostatistical approach for mapping thematic classification accuracy and evaluating the impact of inaccurate spatial data on ecological model predictions
Pöthkow et al. Nonparametric models for uncertainty visualization
Ma et al. An efficient two-stage sampling method for uncertainty quantification in history matching geological models
Fernandez Martinez et al. Reservoir characterization and inversion uncertainty via a family of particle swarm optimizers
Jafarpour et al. Transform-domain sparsity regularization for inverse problems in geosciences
Tompkins et al. Scalable uncertainty estimation for nonlinear inverse problems using parameter reduction, constraint mapping, and geometric sampling: Marine controlled-source electromagnetic examples
Lee et al. Bayesian inversion with total variation prior for discrete geologic structure identification
Arnold et al. Optimisation of decision making under uncertainty throughout field lifetime: A fractured reservoir example
Fernández-Martínez Model reduction and uncertainty analysis in inverse problems
Lindberg et al. Blind categorical deconvolution in two-level hidden Markov models
Rimstad et al. Approximate posterior distributions for convolutional two-level hidden Markov models
WO2011159255A2 (fr) Analyse de données à grand nombre de dimensions
Zhang et al. An introduction to variational inference in geophysical inverse problems
Xiao et al. Surrogate-assisted inversion for large-scale history matching: Comparative study between projection-based reduced-order modeling and deep neural network
Meles et al. Bayesian tomography with prior-knowledge-based parametrization and surrogate modelling
Mahjour et al. Evaluation of unsupervised machine learning frameworks to select representative geological realizations for uncertainty quantification
You et al. Computing the finite time Lyapunov exponent for flows with uncertainties
Friedel et al. Hybrid modeling of spatial continuity for application to numerical inverse problems
Zhang et al. Bayesian spatial modelling for high dimensional seismic inverse problems
Pradhan et al. Consistency and prior falsification of training data in seismic deep learning: Application to offshore deltaic reservoir characterization
Kang et al. A hierarchical model calibration approach with multiscale spectral-domain parameterization: application to a structurally complex fractured reservoir

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 11728436

Country of ref document: EP

Kind code of ref document: A2

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 11728436

Country of ref document: EP

Kind code of ref document: A2