WO2011160201A1

WO2011160201A1 - Method and system for determining an estimated model parameter from a measured field

Info

Publication number: WO2011160201A1
Application number: PCT/CA2011/000707
Authority: WO
Inventors: Felix Hermann; Eldad Haber
Original assignee: University of British Columbia
Current assignee: University of British Columbia
Priority date: 2010-06-21
Filing date: 2011-06-15
Publication date: 2011-12-29
Anticipated expiration: 2012-12-21

Abstract

The present disclosure is directed at a method and system for determining an estimated model parameter from a measured field. The measured field may be generated using, for example, multiple sources located through a domain of interest, such as a certain volume of land, and measured using, for example, multiple receivers also located within the domain of interest. The method includes casting an inverse problem relating the measured field to the model parameter as a stochastic optimization problem, in which the stochastic optimization problem includes minimizing an objective function having a misfit expressed as an expected value of a function that depends on a random variable; iteratively optimizing the objective function using a stochastic optimization technique; determining when the objective function is reduced to a desired value; and when the objective function is reduced to the desired value, identifying a model parameter that corresponds to the objective function at the desired value as the estimated model parameter. Exemplary model parameters that can be estimated include conductivity and seismic velocity. Beneficially, using stochastic as opposed to deterministic techniques can result in savings in computational resources in the forms of FLOPS and memory.

Description

METHOD AND SYSTEM FOR DETERMINING AN ESTIMATED MODEL

PARAMETER FROM A MEASURED FIELD

TECHNICAL FIELD

[0001] The present disclosure is directed at a method and system for determining an estimated model parameter, such as a conductivity model, from a measured field, such as a static electric field. More particularly, the present disclosure is directed at a method and system for determining an estimated model parameter by employing stochastic optimization instead of deterministic optimization, resulting in a technical benefit of savings in computational resources.

BACKGROUND

[0002] Solving parameter estimation problems that employ partial differential equations as constraints and that have multiple right hand sides is performed relatively commonly in applications such as those related to impedance tomography, DC resistivity, electromagnetic imaging, seismic imaging, and hydrology. Such problems are typically mathematically expressed as follows: inin J(yi , . . . , y_n , u) ² + R(u)

s.t ¾(¾, «) = A{u)y_j - ¾· = 0 j = 1 , . . . , N where u e R^m is the discretized model (or control) and y_j e R" j =1 , ... , N is the discretized field (or state). The constraints C_j(v_j,u) = 0 are discretized partial differential equations that share the same model u and are linear with respect to qj. A(u) is presumed to be invertible for all relevant u. The objective function J is composed of two functions: the first is the misfit that measures the fit of y,, projected by the matrix P to some given data dj; and the second, R(u), is assumed to be a convex regularization functional. Every row in P represents a measurement at a certain point in space and time.

[0003] Solving Equation (1) conventionally is very computationally intensive in terms of both FLOPS and memory. Accordingly, research and development continues into methods and apparatuses for solving such parameter estimation problems in more computationally efficient ways.

SUMMARY

[0004] According to a first aspect, there is provided a method for determining an estimated model parameter from a measured field. The method includes casting an casting an inverse problem relating the measured field to the model parameter as a stochastic optimization problem, wherein the stochastic optimization problem comprises minimizing an objective function comprising a misfit expressed as an expected value of a function that depends on a random variable; iteratively optimizing the objective function using a stochastic optimization technique; determining when the objective function is reduced to a desired value; and when the objective function is reduced to the desired value, identifying a model parameter that corresponds to the objective function at the desired value as the estimated model parameter.

[0005] The measured field may intersect a medium. The method may accordingly also include reconstructing an image of the medium from the estimated model parameter.

[0006] The method may also include generating the measured field using a source, and measuring the measured field using a receiver. Different sources and receivers may be used in different aspects. For example, to perform seismic tomography, the sources may be placed in boreholes and generate a pressure field or a displacement field in response to an acoustic signal, the receivers may be placed on the surface, and the estimated model parameter may be seismic velocity. In another aspect and potential application, the sources may generate a magnetic field in response to a current signal and the estimated model parameter may be conductivity. In another aspect and potential application, the sources may generate a static electric field in response to a current signal, and the estimated model parameter may be conductivity. Different applications are possible in a variety of aspects so long as estimating the model parameter from the measured field is an inverse problem that can be cast as a stochastic optimization problem according to the aspects described herein.

[0007] At least one of the receivers may not measure data from at least one of the sources. This may occur because of relative motion between the sources and receivers. [0008] The stochastic optimization technique applied to the stochastic optimization problem may be sample average approximation, and alternatively may be stochastic approximation.

[0009] According to another aspect, there is provided a system for determining an estimated model parameter from a measured field. The system includes a processor; and a memory communicatively coupled to the processor and having encoded thereon statements and instructions to cause the processor to execute a method, which includes casting an inverse problem relating the measured field to the model parameter as a stochastic optimization problem, wherein the stochastic optimization problem comprises minimizing an objective function comprising a misfit expressed as an expected value of a function that depends on a random variable; iteratively optimizing the objective function using a stochastic optimization technique; determining when the objective function is reduced to a desired value; and when the objective function is reduced to the desired value, identifying a model parameter that corresponds to the objective function at the desired value as the estimated model parameter.

[0010] The system may also include a display communicatively coupled to the processor, and the method may also include reconstructing an image of a medium that the measured field intersects from the estimated model parameter.

[001 1] The system can also include one or more sources, communicatively coupled to the processor and located at source positions, configured to emit the measured field when stimulated using a signal; and one or more receivers, communicatively coupled to the processor and located at receiver positions separated from the source positions, configured to measure the measured field. Different types of sources and receivers can be used. For example, in one aspect and to perform seismic tomography, the sources may be placed in boreholes and configured to generate a pressure field or a displacement field in response to an acoustic signal, the receivers may be placed on the surface, and the estimated model parameter may be seismic velocity. In another aspect and potential application, the sources may be configured to generate a magnetic field in response to a current signal and the estimated model parameter may be conductivity. In another aspect and potential application, the sources may be configured to generate a static electric field in response to a current signal, and the estimated model parameter may be conductivity. Different applications are possible in a variety of aspects so long as estimating the model parameter from the measured field is an inverse problem that can be cast as a stochastic optimization problem according to the aspects described herein.

[0012] At least one of the receivers may not measure data from at least one of the sources. This may be due, for example, to relative motion between the sources and receivers.

[0013] The stochastic optimization technique applied to the stochastic optimization problem may be sample average approximation, and alternatively may be stochastic approximation.

[0014] According to another aspect, there is provided a computer readable medium having encoded thereon statements and instructions to cause a processor to execute a method of any of the foregoing aspects or any useful combination thereof.

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] In the accompanying drawings, which illustrate one or more exemplary embodiments:

[0016] Figure 1 depicts a system for determining an estimated model parameter from a measured field, according to one embodiment.

[0017] Figure 2 depicts a method for determining an estimated model parameter from a measured field, according to another embodiment.

[0018] Figures 3(a) - (d) depict approximations to the misfit using different numbers of vectors, which can form part of the embodiments of Figures 1 and 2.

[0019] Figure 4 depicts a schematic of a system used to perform electromagnetic impedance tomography, according to another embodiment.

[0020] Figures 5(a) - (d) depicts results of the electromagnetic impedance tomography performed using the system of Figure 4. [0021 ] Figure 6 depicts a seismic tomography experiment in which a model parameter to be estimated is seismic velocity, and which shows both sources and receivers.

[0022] Figures 7(a) - (f) depict results of the seismic tomography experiment of Figure 6.

DETAILED DESCRIPTION

[0023] Directional terms such as "top," "bottom," "upwards," "downwards," "vertically" and "laterally" are used in the following description for the purpose of providing relative reference only, and are not intended to suggest any limitations on how any article is to be positioned during use, or to be mounted in an assembly or relative to an environment.

[0024] In several fields, such as impedance tomography, DC resistivity, electromagnetic imaging, seismic imaging, and hydrology, tangible phenomena can be imaged by performing parameter estimation. The parameters that are estimated are those of a model, such as a conductivity model, for a particular domain of interest. Such parameter estimation can be expressed mathematically as Equation (1), which is reproduced below: min J{yi, . . . , y_n, u) =—∑\\ Py_j - d₃ f + R(u)

3

s.t Cj (¾- , u) = A( )y_j - ¾^■ = 0 j = 1 , . . . , N

[0025] For example, parameter estimation may be performed to estimate values for conductivity or susceptibility of a target, such as an oil deposit. The domain of interest may be a certain volume of land. The measured field, y, used to estimate u may be an electromagnetic field. The electromagnetic field may be generated using multiple current-injected sources located within boreholes in the land, and the electromagnetic field may be detected using receivers located on the surface of the land, for example. The electromagnetic field is affected if it intersects an oil deposit. By estimating values for the conductivity/susceptibility model from measured values for the electromagnetic field, it is accordingly possible to determine where the oil deposit is by identifying the appropriate fluctuations in the conductivity/susceptibility model. Similar methods can be employed in other fields, such as in medicine to image bones by estimating elastic parameters, to obtain images of certain media of interest. [0026] Conventionally, Equation ( 1 ) is solved deterministically. "Deterministically" in this context refers to solving Equation (1) using algorithms that behave predictably; i.e., given a particular input, a deterministic algorithm will always produce the same output. In contrast, the embodiments described herein instead solve Equation (1) stochastically. "Stochastically" in this context refers to solving Equation (1) using algorithms that behave, at least to a certain degree, randomly. Solving Equation (I) stochastically as opposed to deterministically can result in a technical benefit of significant savings in computational resources (e.g.: savings in one or both of FLOPS and memory used). Compared to conventional deterministic methods of parameter estimation, the stochastic embodiments described herein can be about ten times faster and also typically consume significantly less memory. This results in savings in any one or more of time, hardware (processing and memory requirements), and money. For example, when performing a single instance of parameter estimation in the oil and gas industry, the embodiments described herein can result in savings of roughly $100,000 to $1,000,000 relative to conventional methods.

[0027] Referring now to Figure 1 , there is depicted a system 100 for determining an estimated model parameter from a measured field; the system 100 can be used, for example, to estimate DC resistivity. The system 100 includes a surveying truck 102 that is stationed over a particular volume of land to be surveyed. One goal of the system 100 is to be able to locate the presence of the medium 1 14 located within the land. Located in boreholes within the land are sources 104 in the form of permanent current electrodes. Also positioned in the land is a receiver 106 in the form of a cone-mounted potential electrode. The cone-mounted potential electrode records voltage differences relative to a potential electrode 105, located on or near the surface of the land. Both the sources 104 and the receiver 106 are communicatively coupled to, and controlled via, the truck 102. The system 100 is accordingly configured to generate and to measure a field in the form of a static electric field. The static electric field intersects with the medium 1 14, which affects the field and which accordingly allows the medium 1 14 to be identified.

[0028] Measurements obtained via the truck 102 are sent to a processor 108 for analysis, where a log conductivity model for the volume of land can be estimated. The processor 108 is communicatively coupled to a memory 1 10, which contains statements and instructions to cause the processor 108 to perform a method for parameter estimation, such as the one described with reference to Figure 2, below. A display 1 12 is also communicatively coupled to the processor 108 to allow the processor 108 to resolve and output an image identifying the location of the medium 1 14.

[0029] In contrast to attempting to solve Equation (1) deterministically, the embodiments described herein solve Equation (1 ) stochastically; Equation ( 1) is cast as a stochastic optimization problem involving estimation of the trace. To arrive at this solution, Equation (1), unconstrained, can first be re-written as: min J(u) = pA{u)-^lQ - Dt| + R(u)

u^{' ' 2} (2) where Q = (qi, qN), D = (di, d_N), and ||-||_F denotes the Frobenius norm. The difficulty is to evaluate the residual term S(u) = PA(U)^"'Q - D and its norm. To address this difficulty, a stochastic interpretation of the trace can be used.

[0030] For a random variable/vector w with 0 mean and covariance matrix I, the following holds:

E„ (w^TS(u) ^TS(u)w) = trace(S(u)^TS(u)) =

[003 1] Using Equation (3) the deterministic problem of Equation (2) can be replaced with a stochastic optimization problem, as follows: min E_w J

[0032] In Equation (4), E_w is the expected value with respect to w. While Equation (4) is mathematically identical to Equation (2), it requires fewer computational resources to solve. Because Equation (4) is in the form of a "stochastic optimization problem", a variety of algorithms referred to herein as "stochastic optimization algorithms" can be used to solve Equation (4). As long as the evaluation of the expectation in Equation (4) can be done efficiently and with sufficient accuracy, there is an efficiency gain over using deterministic algorithms. [0033] In the stochastic formulation, each realization Wj of w involves a single solve for the partial differential equation; i.e.:

(PA(u) 'Q - D)w_j = PAiuV'tQw_j) - Dw_j, which implies that given Wj, only a single partial differential equation is solved in order to evaluate the stochastic misfit function. This is in contrast to the N partial differential equations to be solved when evaluation Equation (2). Accordingly, relatively efficient algorithms and methods can be employed.

[0034] In particular, two end members of stochastic optimization algorithms can be considered: Stochastic Average Approximation (SAA) and Stochastic Optimization (SA).

[0035] When using SAA, one approximates the expected value by a Monte-Carlo approximation using the sum:

E_w J( i, w) i¾— ,/(«, Wj ).

j= i « N different realizations are used to approximate the expected value.

[0036] When using SA, a single realization is used and the gradient V_uJ(u,Wj) is determined. A step in the negative direction of the gradient is then taken, averaged with the previous steps and the process repeats. For each step of the SA embodiment, only a single realization is used. Applying this to solving Equation (4) implies that only a single partial differential equation and adjoint are solved at each iteration, rather than 2N in Equation (2) and 2 when using SAA. Both SA and SAA can lead to increased performance relative to conventional methods because each approximation may require fewer partial differential equation solves. Applying both SA and SAA are considered in more detail, below.

Solution Through SAA

[0037] When using SAA, the expectation can be approximated as follows: min E_w (J_w(u, w)) « J_K{u; w . . (u).

Even with K = 1 , the evaluation of the trace has been found to be sufficiently accurate in some contexts.

[0038] In order to choose the vectors Wj stochastic trace estimators are considered. If w is 0 mean with covariance I then E_ww^TS^TSw = trace(S^TS). The expected value is identical for all w's drawn from any distribution with 0 mean and covariance I. However, as approximations are being performed, the goal is to obtain this mean with a relatively small variance, and in certain embodiments the smallest variance possible; this can help in generating a relatively small, and in certain embodiments the smallest, number of right hand sides in the partial differential equation. This can be obtained when w is a random vector drawn from the distribution of +/- 1. w is accordingly chosen from this distribution. In alternative embodiments, other choices of w are possible, although they are suboptimal in that their variance is higher.

[0039] SAA beneficially separates the questions of approximating the expectation and the optimization algorithm to be used. In order to solve the optimization problem any robust algorithm can be used. The LBFGS method is applied to the unconstrained problem in the embodiments explicitly discussed herein; however, in alternative embodiments, other methods may be applied to solve the problem.

Solution through SA

[0040] An exemplary algorithm in pseudocode that can be used for SA is given below:

1 : Initialize the solution u = uo

2: while not converge do

3: choose a realization Wk

4: approximately solve the optimization problem s = arg mins_u J(uk + 8u, Wk) 5: set Uk+ i = Uk + ys

6: average ι¾+ι = [l/(k+l)](∑jiij + Uk+i)

7: end while

[0041] Considering the foregoing, two questions arise: first, how to pick an appropriate realization Wk; and second, how to approximately solve the optimization problem and choose the step size γ in step 5 of the foregoing algorithm. As discussed above, a beneficial choice for w is obtained by choosing w to be a random vector of {+/-1 } . For an optimization algorithm, a steepest descent type method may be used, and potentially Gauss-Newton and LBFGS methods as well. As above, to facilitate comparison with SAA results, the exemplary embodiments discussed herein utilize the LBFGS method.

The Quality of the Approximation

[0042] The quality of the approximation discussed above depends on the stochastic approximation of Equation (3). For statistical problems which are set appropriately, convergence is guaranteed as the number of samples goes to infinity. Since Monte Carlo methods are used herein, convergence is typically related to the square root of the number of samples. The error is therefore associated with the fact that the problem is approximated using a finite number of sampling points and the variance of the solution is accordingly the relevant quantity to examine. To illustrate this principle, an exemplary DC resistivity problem is used. In such a DC resistivity problem, the matrix A results from discretization of the problem:

V · exp(u) Vy_J = q_j V¾ - n = 0 / ;¾ = 0 j = l, . . . , N where the goal is to recover log conductivity u given measurements of the potential fields y_j, which result from the sources 104 qj. Modern acquisition systems can have hundreds and thousands of the sources 104 spread over a physical grid and in boreholes. The partial differential equation shown in Equation ( 1) results from performing a finite volume or finite element discretization of the differential equation. [0043] The goal is to approximate the misfit(u) = PA(u)-^lQ - D\\% using stochastic approximation. To begin, u is set to be constant, e.g. 10^"2, the gradient of the misfit, s = V_umisfit, is computed, and a ID function is defined as follows: f( ) := misfit ( it + as)— \\ PA(u + as) ¹Q— D

[0044] The following approximation is then used: f_K (a) = misfit(u + as; wi , w_K)—— || (P_4(u -I- as) ¹Q - D)w_j \\².

[0045] Presuming 500 of the sources 104, the partial differential equation is discretized using 32³ cells, where y and u are discretized on cell-centers. Thus, Q : R⁵⁰⁰ - R^{32 3} and to compute the misfit in the original formulation, A(u) is inverted 500 times. Random realizations [W], ... w_k], where K = 1 , 10, 25, 50 are then experimented with. Evaluating the misfit in this way is roughly a factor of 500 ^"1 cheaper when comparing floating point operations compared to the original evaluation. Figures 3(a) - (d) are plots of the estimated misfit misfit as a function of a. As each approximation depends on a particular choice of w that is random, five different realizations of w are used for each approximation. Although the functions are not identical, when using a single vector to evaluate the misfit, the minimum with respect to a is not far from the minimum obtained by exact evaluation of the misfit.

A Quantitative Evaluation

[0046] In the following discussion, W[i _K] refers to a sample of K vectors [wi,...,w_K].

W'[I K] means the 1^th sample of vectors, u* is the minimizer of the approximate objective function. Accordingly, for a given W^KJ: u* = argrnin J (it: u_(1>A1) =— Y j(u;

J

[0047] Let u^† be the minimum of the exact objective function:

«^† = argmin J(u) = E_w(J_w(u, w))

u and

J(u^f) < J(u*).

[0048] The quality of the approximation can be measured by the gap: gap = J(0 - J(M^†) which can be estimated statistically.

[0049] To estimate J(u*), w jq is resampled and the value of J(U*,W[I ]) is computed for different realizations of Wfi q- If W[i _K] is sampled L times, the mean and the variance can be estimated, obtaining:

J(u*) « j∑ J_K (u*:wf _hK]) = J_L*.

t

An optimization problem does not need to be solved to compute this estimate.

[0050] u can then be fixed to equal u* and the samples W[i _] can be changed, and given the realizations the variance of J(u*,W[i «_.]) can also be estimated:

[0051] The combination of the mean and the variance provides a lower and an upper estimate on J(u*). In particular, the 66% confidence interval is given by:

Jt - σ^* < J(u*) < t + σ ₍₆₎

[0052] More computational work is performed to evaluate J(u^†). Consider, as before, independent and identically distributed samples w'_{[l iK]}, 1 = 1,...,T, and solve T optimization problems: iif = ar g mm JK ( u ; w f_t i ) , = 1 , . . . , T

u

[0053] T different solutions [ui^*,...,u_T*] are obtained. The variance of the function J(u^†) is determined as opposed to its exact value. Using the solutions [UI*, . . .,UT*] the mean and variance of J^† can be determined. First, let

J be the mean of the function value at the different solutions, and let

be the variance of the function of the different solutions. The estimated 66% confidence interval is then given by

4 - < ^j(_M ^f) < + ₍₇₎

[0054] Combining Equations (6) and (7) allows an estimate of the quality of the solution to be obtained. Unlike evaluating J(u*) where only function evaluations are performed, evaluating J(u^†) involves solving some optimization problems. However, in many cases very few optimization problems are solved (e.g.: between 5 and 10), and to efficiently use computational resources these optimization problems are independent and can be solved in parallel, where a suitable system is available. Additionally, from a computational point of view solving these optimization problems can be done quickly if the solution u* of the first optimization problem is a starting point.

Embodiment of Different Observation Operators

[0055] The foregoing embodiments assume that different experiments represented by <¾ have the same observation matrix P. Although many problems share this property this is not always the case. Accordingly, the following embodiments are directed at situations in which P differs for different experiments. Practically, this may represent a situation in which the sources 104 and the receivers 106 move while measurements are taken (e.g.: if the receivers 106 are dragged by a boat during measuring). More generally, this occurs when at least one of the receivers 106 measures no data corresponding to at least one of the sources 104. In these embodiments, the misfit is first rewritten in general as misfit(tt) = \\C Θ (PA(ur^lQ - D)]|| - \\C ø S(u) \\

(8) where as before S(u) = PA(u)^"'Q - D is the residual matrix. Here Θ is the point-wise Hadamard product. The matrix C is, in general, a dense matrix. To illustrate that the foregoing is indeed a generalization, consider the following cases:

• If C = 1, that is, populated by I s, then the previous case results in which all the sources 104 share all the receivers 106.

• If C is dense and populated by different positive numbers, typically each datum has a different standard deviation and Cy is the inverse of the standard deviation of each datum.

If C is sparse then only a subset of source-receiver configuration is chosen. • If C is diagonal, which is a relatively extreme case, it is implied that each of the sources 104 has a single one of the receivers 106. This configuration is relatively common in applied geophysical surveys in which the sources 104 and the receivers 106 move simultaneously.

[0056] As before, the misfit is represented by the following expectation: misfit(u) =

= E_w ø S(tt)M|^a) = E» (II (C O (ΡΑ(^η)^~^ - D)) wf) .

[0057] Evaluating Equation (9) can be difficult. Recall that for C = 1:

(1 Θ S(u)) w = S(u)w = PA(u)^~l(Qw) - Dw

[0058] However, when C does not equal 1, the same methodology cannot be followed.

The difficulty is that the Hadamard product does not associate with matrix-vector product and therefore one cannot compute the product S(u)w first and only then compute the Hadamard. This provides a challenge in extending the method discussed in respect of the embodiments above, in which P is constant for all different experiments. The following discusses two different approaches for solving this problem.

[0059] In some applications it is possible to decompose C into low rank matrices, or at least, to approximate C by a product small rank matrices. Consider the case where

C = YZ^T ( 10) where Y and Z are n_r x k and n_s x k matrices respectively with k « min (n_r,n_s). Using the low rank representation of C the following can be proved.

[0060] First, let R be an n_r x n_s matrix with columns R, i = l,...,n_s and let Y and Z be n_r x k and n_s x k matrices respectively, with columns Y_j and Z_j, j = l,...k. Finally, let w be a vector of size n_s. Then: n_r k k

(C Θ R)w = (YZ^T Θ R)w =∑ J¼ 0∑ ^Yi^z» =∑∑ ^M¾¾ © 31

i= l 3 j = l i-l

The foregoing implies that it is possible to obtain a stochastic representation to the

J] ¾ Θ (P^TA(u ¹Q - D)(Zi © w)

(1 1)

[0062] The above implies that it is possible to obtain a stochastic approximation similar to the case that C = ee^T using k matrix solves where k is the rank of C. In cases that C is not low rank, it is possible to approximate C by a low rank matrix and use the above decomposition as a cheap alternative to the original misfit function.

The General Case

[0063] A number of different ways exist to obtain a stochastic approximation to the trace in this case. Two such approximations are discussed in the following.

[0064] One way is based on direct estimation of the product C Θ S(u). It is straightforward to observe that if w is a random variable with 0 mean and identity covariance then

C Θ S(u) = Ew (diag (S(u)w) C diag (w)) (12) and therefore one arrives to the following stochastic representation of the misfit: misfit(tt) = \\C Θ S(u) _F = \\E_W (diag (S(u)ti/) C diag («/)) I (13) and its approximation:

misfit(' )

[0065] A difference between the approximation of Equation (14) and those discussed in respect of other embodiments discussed above is that in Equation (14), the expected value is inside the norm, and hence the expectation (approximated by the sum) is first evaluated inside and only then is the norm evaluated. Again, only the matrix vector product S(u)w_j is needed in order to evaluate the misfit. Thus the number of partial differential equations to be solved is equivalent to the number of stochastic realizations.

[0066] To obtain the derivative of the approximation to the misfit, vec (diag (S^w^ C diag («¾)) = C_w.S^w where the matrix Cwj is defined by

C_w . = [biockdiag((Cdiag (u¾-))i-, (Cdiag (w_J))_2→, . . . , (Cdiag («¾))*_, )] ^Γ , where the notation 1 - implies the 1^th row of the matrix. Using this notation the approximated misfit can be rewritten as

misfit{w) C_wj iPM i) ⁱQw_j - Dw_j)

[0067] This representation allows for simple differentiation of the approximated misfit

V_ttmisfit(«)

where the matrix G_j = G(Qw_j,u) is defined by G_j(y_j,u) = Vu[A(u)yj]. This formulation of the problem requires K solves of the adjoint problem.

[0068] Other stochastic approximations to the trace can be obtained. Starting from

Equation (9), the following identity can be applied:

(C Θ S( )) w = diag (C diag (w) S(u)^T) and therefore misfit(u) = \\C 0 S{u) f_F = E_w (||diag (C diag (w) S(u)^J)f) ( 15)

[0069] Identity ( 15) does not solve the problem of computing either A(U)^"'Q or A(u)^~TP.

The point here is that evaluating the diagonal of the matrix product

Cdiag (ti/) S(u)^T = C diag («/) {ΡΑ{^η)^~^ - D) ^T is needed without computing the matrix itself. Thus, another stochastic process is used. If x is a random variable for zero mean and covariance Cov(x) = I, then: diag (C diag (w) S{u)^T) = _x (x ø (C diag (w) S(u)^T)x) ,

(16) [0070] Thus the misfit can be replaced with a double stochastic process misfit(_u) = i|C ø = E_w (|J¾ (x © (C diag H S(u)^T)x) f) . _{( l y)}

[0071] Again, using Monte-Carlo to obtain an approximation results in:

which requires the computation of S(u)^Tx that implies the solution of the adjoint problem. Experimentally, it has been found that in general, the approximation (18) requires more realizations than the approximation (14) in order to achieve similar accuracy of the misfit.

[0072] Using the stochastic misfit approximation misftt(u) given by either (14) or (18) the SAA method can be used to solve the stochastic optimization problem and approximately recover the model u. Notably, since the expected value is inside the norm using SA is not justified.

Numerical Experiments

[0073] The following discusses numerical experiments performed to test the foregoing approximations and to experiment with different methods and parameters. Both the exemplary DC resistivity problem, discussed above, is used, as well as a seismic tomography/full-wave inversion problem. In the seismic tomography problem, consider the case where the matrix A results from discretization of the problem

(Δ + fc²(ti)) if, = <¾, in Ω C E², ¾- = 0 on Γ = 0Q_} j = 1, . . . , N. where k²(u) = 2nfu is the wave number with f the temporal frequency, and u = c^"2 with c the acoustic wave-speed. In this case the goal is to recover the velocity c = u^"1 given a measurements of the wave-field fields yj for each source experiment. Modern acquisition systems can have hundreds and thousands of sources spread over a physical grid and in boreholes.

Experiments with DC Resistivity Problem

[0074] For the DC resistivity problem, the code for the forward problem and Jacobians is based on a software package described at A. Pidlisecky, E. Haber, and R. Knight, Resinvm3d: A matlab 3-d resistivity inversion package, Geophysics, 72:212-224, 2006.

[0075] In particular, the following compares conventional deterministic methods with the stochastic methods of the embodiments discussed herein. Since the deterministic methods utilize large amounts of computer memory, the LBFGS algorithm is used for all methods such that the comparison between deterministic and stochastic methods is on equal footing.

[0076] The following utilizes a mesh of 32³ cells. For the embodiment in which all the sources 104 share the same receivers 106, the experiments are performed assuming 250 each of the sources 104 and the receivers 106. As depicted in Figure 4, which is a schematic of the system used to perform the DC resistivity experiment, the sources 104 are located in boreholes while the receivers 106 are located in a grid-like arrangement on the surface. Also as depicted in Figure 4, the media 1 14 are an ellipsoid conductor with a conductivity of 10 S/m and a resistor with a conductivity of 0.1 S/m. For simplicity a smooth recovery is used, and

R(u) = £u^TVl _ftti, where h is a discretization of the gradient operator and p is chosen such that a fit to the data is obtained based on the discrepancy principle.

Recovery When All the Sources 104 Share the Same Receivers 106

[0077] In a first experiment in which it is assumed that all the sources 104 share all the same receivers 106, the stochastic problem of Equation (4) results. To avoid the "inverse crime", data for the true model is generated using a fine grid of 64³. Three different methods are then used to solve the inverse problem; i.e., to solve for the model u based on the measured field. In the first method, the deterministic LBFGS method is used to yield what is conventionally considered an "exact" solution of the stochastic optimization problem. In a second experiment and according to one of the embodiments, five vectors w are selected and the SAA embodiment is used to stochastically obtain a solution. In a third embodiment, the SA embodiment is used with a single realization at each step, while replacing realizations as the method progresses. Figures 5(a) - (d) show the different results: Figure 5(a) shows the true solution; Figure 5(b) shows the deterministic solution; Figure 5(c) shows the solution using SA; and Figure 5(d) shows the solution using SAA. The results obtained using the stochastic optimization methods according to the depicted embodiments closely resemble the deterministic optimization: the relative difference in the recovered models by the different techniques is below 3%. In this particular example, it is possible to compute the true value of the misfit and its stochastic approximation. The relative difference between the misfit obtained by the deterministic optimization method and the misfit obtained by the SAA embodiment was 0.1 1%, and the difference between the misfit obtained by the SA embodiment and the misfit obtained by the deterministic optimization was 0.13%. Thus, the SA and SAA embodiments are able to closely approximate the true solution, while being much more computationally efficient in terms of computational resources used, as illustrated below in Table 1 by the lower number of partial differential equations solved while performing the SA and SAA embodiments:

Table 1: Number of Iterations and Partial Differential Equation Solves Used by

Deterministic, SAA, and SA Embodiments

[0078] Table 1 shows that a factor of about 60 in terms of partial differential equation solves is observed between the stochastic and deterministic methods. While the estimates obtained using the different methods are roughly similar, the savings in resources provided by the SA and SAA embodiments allows larger problems to be tackled using a fraction of the computational effort used when solving deterministically.

[0079] For the foregoing:

• For the deterministic and SAA optimization, computations are stopped when the norm of the gradient is 10^"3 of its initial value. • For the SA embodiment, LBFGS is used, and stopped when the change in the solution is smaller than 10^"3.

• For the SAA embodiment, five vectors are used to obtain the solution. The estimates discussed above are used as stopping criteria. The number of w's is determined using experimentation. When evaluating the gap = J(u^*) - J(u^†) the problem is first solved for five different samples starting from the solution obtained by the first sample. For all problems the solutions for different samples converge in a single iteration, which indicates that iterations began close to the solution and the variability between solutions is small.

The question of how many realizations of w's to use is determined through experimentation. From the starting point, the variability of the misfit is experimented with. For a small set of w's (for example using a single vector), a reasonable approximation to the misfit (and a reasonable recovery of u) is obtained, but repeating the experiment a non-negligible variance of about 7% between the solutions is found from different realizations. When using five realizations the results change relatively little (less then 3%) between different realizations. Determining the number of realizations a priori can be beneficial.

Recovery When All the Sources 104 Do Not Share the Same Receivers 106

[0080] In a second set of experiments, the stochastic approximations given in Equations

( 14) and (18) are used to solve the problem where the matrix C is not equal to 1. Three different cases are tested:

• First, the matrix C is dense and contains the standard deviation of each datum.

1 % noise is assumed, and thus each entry in C is roughly 1 % of the value of the corresponding datum.

• Second, the matrix C is sparse, made of 0's and l 's, and contains only 10% nonzero entries. The entries are randomly chosen. • Third, there are an equal number of the sources 104 and the receivers 106 and the matrix C is diagonal, C = I; i.e., each of the sources 104 has a single one of the receivers 106. This case is relatively common in certain geophysical applications in which the sources 104 and receivers 106 simultaneously move.

[0081] Three different optimization problems are solved for each of the foregoing:

• one, a deterministic LBFGS;

• two, SAA with the stochastic approximation in Equation ( 14) (referred to as SAA1 in Table 2, below); and

• three, SAA with the stochastic approximation in Equation (18) (referred to as SAA2 in Table 2, below).

[0082] In this second set of experiments, 500 of the sources 104 and 500 of the receivers

106 are used in the same configuration as shown in Figure 4. Again the number of iterations and partial differential equation solves are recorded for each of the cases, as well as the difference between the models u recovered by stochastic optimizations relative to the models recovered for the deterministic one. The results for the different cases are summarized in Table 2, below.

Table 2: Number of Iterations and Partial Differential Equation Solves Used by the Dense,

Sparse, and Diagonal Cases

SAA1 32 325 1.7%

SAA2 36 745 5.3%

Diagonal Deterministic 27 28,500 0%

SAA1 32 331 0.9%

SAA2 39 789 6.2%

[0083] In respect of the experiments generating the data for Table 2:

• For SAA1, five vectors were used and experimentally it was found that this was sufficient to obtain results with relatively little variability.

• For SAA2, two vectors were used for wi and five vectors for w₂, and the computations are almost double compared to SAA1. As can be seen in Table 2, SAA1 results in a smaller variance compared to SAA2.

• All the optimization algorithms were terminated when the gradient was reduced to 10^"3 of its original value.

[0084] As shown in Table 2, the model u recovered by statistical embodiments is similar to the models recovered by the deterministic models. Also as above, a technical benefit is that the depicted embodiments are performed in a fraction of the computational cost (e.g.: measured in FLOPS) of the deterministic methods.

Experiments with Seismic Tomography

[0085] In this experiment, the 2D seismic borehole tomography example presented above is employed. Two boreholes are assumed to be 1 meter apart and it is assumed that the left borehole contains 500 of the sources 104 while the right borehole contains 500 of the receivers 106. A frequency of 2 Hz is used that, for the material properties assumed in this experiment, implies roughly 3 wavelengths between boreholes. A sketch of the seismic tomography experiment and the true model is plotted in Figure 6, in which the estimated model parameter is seismic velocity. The sources 104 are represented by the dark bar on the right of Figure 6, while the receivers 106 are represented by the light bar on the left of Figure 6.

[0086] In this experiment, SAA is used; however, unlike in the DC resistivity experiment where the number of vectors is set a priori, continuation is used to determine the number of samples needed. A single sample is used to start, the optimization problem is solved, and then two new samples are used starting from the previous solution. This is continued until the recovered model does not change much (e.g.: less than 1% change). Figures 7(a) - (f) plot the models obtained following this process using one to six samples.

[0087] As can be seen, using a single sample was not particularly beneficial; however, when using five and then six the result was relatively unchanged.

[0088] In the foregoing embodiments, although LBFGS is used, other algorithms can be used that converge relatively better and faster. For example, the Gauss-Newton method and methods based on equality constrained optimization can potentially be used and should have better numerical properties. Using better optimization techniques may further expand the gap between stochastic and deterministic optimization techniques.

[0089] Given the foregoing, and referring back to the system 100 of Figure 1, the following is a discussion of how the processor 108 handles the measurements it obtains using the receivers 106. The processor 108 performs a method 200, depicted in Figure 2, for determining an estimated model parameter from a measured field, according to another embodiment. The method beings at block 202, from which the processor 108 transitions to block 204 where it casts the inverse problem (i.e.: a full-waveform inversion problem) as a stochastic optimization problem so as to be able to apply a stochastic optimization technique, such as SA and SAA, to it as opposed to conventional deterministic algorithms. As discussed above, and as referenced in block 204, employing stochastic optimization involves minimizing an objective function that is present when the inverse problem is cast as a stochastic optimization problem. The processor 108 then proceeds to block 206, where the objective function is iteratively optimized using a stochastic optimization technique such as SA or SAA. Eventually, the objective function will be reduced to a desired value; the desired value may be user controlled, depending on the level of computational resources available and the accuracy of the estimated model parameter desired. When the processor 108 has determined that the objective function has been reduced to the desired value, the processor 108 proceeds to block 210 where the model parameter that corresponds to the objective function reduced to the desired value is adopted, or identified as, the estimated model parameter. The estimated model parameters can optionally be used to generate an image of one of the media 1 14 located in the domain of interest. Once the estimated model parameter is obtained, the processor 108 proceeds to block 212, where the method 200 ends.

[0090] Any of the foregoing methods may be automated using any suitable type of processor, which may include a programmable logic controller, microprocessor, microcontroller, application specific integrated circuit, field programmable gate array, multi-core processor, an assembly of processors configured to execute code in parallel, or the like. Any of the foregoing methods may also be encoded on to a memory or other computer readable medium, which may be communicatively coupled to the processor. The memory may, for example, be any suitable type of semiconductor or disc based memory, such as RAM (whether volatile or non-volatile), ROM, hard disk drives, CD-ROMs, and DVD-ROMs.

[0091] While particular embodiments have been described in the foregoing, it is to be understood that other embodiments are possible and are intended to be included herein. It will be clear to any person skilled in the art that modifications of and adjustments to the foregoing embodiments, not shown, are possible.

Claims

A method for determining an estimated model parameter from a measured field, the method comprising:

(a) casting an inverse problem relating the measured field to the model parameter as a stochastic optimization problem, wherein the stochastic optimization problem comprises minimizing an objective function comprising a misfit expressed as an expected value of a function that depends on a random variable;

(b) iteratively optimizing the objective function using a stochastic optimization technique;

(c) determining when the objective function is reduced to a desired value; and

(d) when the objective function is reduced to the desired value, identifying a model parameter that corresponds to the objective function at the desired value as the estimated model parameter.

A method as claimed in claim 1 wherein the measured field intersects a medium, and further comprising reconstructing an image of the medium from the estimated model parameter.

A method as claimed in any one of claims 1 and 2 further comprising generating the measured field using a source, and measuring the measured field using a receiver.

A method as claimed in claim 3 wherein the source comprises multiple sources and the receiver comprises multiple receivers, and wherein at least one of the receivers does not measure data from at least one of the sources due to relative motion between the at least one of the receivers and the at least one of the sources.

A method as claimed in any one of claims 1 to 4 wherein the stochastic optimization technique comprises sample average approximation.

6. A method as claimed in any one of claims 1 to 4 wherein the stochastic optimization technique comprises stochastic approximation.

7. A method as claimed in any one of claims 1 to 6 wherein the measured field is generated in response to a current signal, the field is a static electric field, and the model parameter is conductivity.

8. A method as claimed in any one of claims 1 to 6 wherein the measured field is generated in response to either a current signal or a magnetic field, the field is a magnetic field, and the model parameter is conductivity.

9. A method as claimed in any one of claims 1 to 6 wherein the measured field is generated in response to an acoustic signal, the field is either a pressure field or a displacement field, and the model parameter is seismic velocity.

10. A system for determining an estimated model parameter from a measured field, the system comprising:

(a) a processor; and

(b) a memory communicatively coupled to the processor and having encoded thereon statements and instructions to cause the processor to execute a method comprising:

(i) casting an inverse problem relating the measured field to the model parameter as a stochastic optimization problem, wherein the stochastic optimization problem comprises minimizing an objective function comprising a misfit expressed as an expected value of a function that depends on a random variable;

(ii) iteratively optimizing the objective function using a stochastic optimization technique;

(iii) determining when the objective function is reduced to a desired value; and (iv) when the objective function is reduced to the desired value, identifying a model parameter that corresponds to the objective function at the desired value as the estimated model parameter.

11. A system as claimed in claim 10 further comprising a display communicatively coupled to the processor, and wherein the method further comprises reconstructing an image of a medium that the measured field intersects from the estimated model parameter.

12. A system as claimed in any one of claims 10 and 1 1 further comprising:

(a) a source, communicatively coupled to the processor and located at a source position, configured to emit the measured field when stimulated using a signal; and

(b) a receiver, communicatively coupled to the processor and located at a receiver position separated from the source position, configured to measure the measured field.

13. A system as claimed in claim 12 wherein the source comprises multiple sources and the receiver comprises multiple receivers, and wherein at least one of the receivers does not measure data from at least one of the sources due to relative motion between the at least one of the receivers and the at least one of the sources.

14. A system as claimed in any one of claims 10 to 13 wherein the stochastic optimization technique comprises sample average approximation.

15. A system as claimed in any one of claims 10 to 13 wherein the stochastic optimization technique comprises stochastic approximation.

16. A system as claimed in any one of claims 10 to 15 wherein the measured field wherein the measured field is generated in response to a current signal, the field is a static electric field, and the model parameter is conductivity.

17. A system as claimed in any one of claims 10 to 15 wherein the measured field is generated in response to either a current signal or a magnetic field, the field is a magnetic field, and the model parameter is conductivity.

18. A system as claimed in any one of claims 10 to 15 wherein the measured field is generated in response to an acoustic signal, the field is either a pressure field or a displacement field, and the model parameter is seismic velocity.

19. A computer readable medium having encoded thereon statements and instructions to cause a processor to execute a method as claimed in any one of claims 1 to 9.