A PENALTY METHOD FOR PDE-CONSTRAINED OPTIMIZATION

• US 20160070023A1
• Filed: 04/23/2014
• Published: 03/10/2016
• Est. Priority Date: 04/24/2013
• Status: Abandoned Application

• Current Assignee: University of British Columbia
• Original Assignee: University of British Columbia
• Inventors: Herrmann, Felix J., Van Leeuwen, Tristan
• Application Number:
  US14/786,146
• Publication Number:
  US 20160070023A1
• Time in Patent Office: Days
• US Class Current: 1/1
• CPC Class Codes: G01V 20/00 Geomodelling in general
• CPC Class Codes: G06F 17/13 Differential equations using digital differential analysers G06F7/64

Abstract

The invention is directed to a computer-implemented method for obtaining a physical model having physical model parameters wherein solutions to one or more partial-differential-equations (PDE'"'"'s) are calculated and wherein (i) an appropriate initial model is selected, (ii) setup a system of equations referred to as the data-augmented PDE for the field, comprising of the discretized PDE, the sampling operator, the source function and the observed data, and solve the data-augmented PDE in a suitable manner to obtain a field that both satisfies the PDE and fits the data to some degree, (iii) setup a system of equations by using the PDE, the source function and the field obtained in step (ii) and solve this system of equations in a suitable manner to obtain an update of the physical model parameters and repeat steps (ii)-(iii) until a predetermined stopping criterion is met.

Claims

1-21. -21. (canceled)

22. A computer-implemented method for obtaining a physical model having physical model parameters, wherein solutions to one or more partial-differential-equations (PDEs) are calculated, wherein the PDEs are constrained by observed data from a system whose response to a known source function and known sampling operator is modeled by solutions of those PDEs for a given physical model using physical model parameters, wherein the method comprises:

(i) selecting an appropriate initial model;

(ii) setting up a first system of equations, referred to as data-augmented PDEs, for a field comprising a discretized PDE, the sampling operator, the source function and the observed data;

(iii) solving the data-augmented PDEs in a suitable manner to obtain a field that both satisfies the one or more PDEs and fits the data to some degree;

(iv) setting up a second system of equations by using the discretized PDE, the source function, and the field obtained in step (ii);

(v) solving the second system of equations in a suitable manner to obtain an update of the physical model parameters; and

(vi) repeating steps (ii)-(v) until a predetermined stopping criterion is met.

23. A method according to claim 22, wherein steps (ii), (iii), and (iv) are performed by minimization of Equation (4) with respect to the fields u_i by solving for each i the data-augmented PDE as defined in Equation (5),

24. A method according to claim 23, wherein steps (iv) and (v) are performed by calculation of a gradient from a modified objective, defined in Equation (6),

φ

_λ

(m)=φ

_λ

(m,ū

(m)), 

 

(6)

wherein ū

is the field obtained in step (ii) and using the expression in Equation (7)

∇

_m φ

_λ

(m)=∇

_mφ

_λ

(m,ū

)=Σ

_i=1^Mλ

²G_i(m,ū

)*(A_i(m)ū

−

q_i). 

 

(7)

and calculation of a sparse approximation of the Hessian (H) by ignoring the (λ

⁴) terms in Equation (8)

∇

² φ

_λ

(m)=Σ

_i=1^Mλ

²G_i*G_i−

λ

⁴G_i*A_i(P_i^TP_i+λ

²A_i*A_i)<sup>−

1</sup>A_iG_i, 

 

(8)

and by summing over all sources, and

calculation of an update for the physical model parameters using the calculated gradient and the calculated inversion of the sparse approximation of the Hessian.

25. A method according to claim 23, wherein the trade-off parameter λ

is increased at a predefined rate ensuring convergence of the algorithm to a solution of the PDE-constrained optimization problem according to equation (1);

min_m,uΣ

_i=1^M1/2∥

P_iu_i−

d_i∥

2/2s.t. A_i(m)u_i=q_i,i=1 . . . M.

26. A method according to claim 23, wherein calculation of solutions of the data-augmented system defined in Equation (5) based on a time-dependent PDE is performed using a time-stepping algorithm.

27. A method according to claim 24, wherein the update for the physical model parameters in step (v) is calculated by solving the data-augmented PDEs of step (iii) in parallel and aggregating the gradient and Hessian across all sources.

28. A method according to claim 24, wherein the number of sources of the source field and/or receivers of the detection operator P, are reduced.

29. A method according to claim 28, wherein the method of reducing the sources comprises a source/receiver-subset selection procedure given by Equation (10)

q_i=Σ

_j=1^Mw_ijq_i, d_i=Σ

_j=1^Mw_ijd_i, 

 

(10)

where i=[1, 2, . . . , M′

] (M′

<

<

M) and w_ij are given by

w_ij=δ

_ij for a randomly chosen jε

[1, M], or by

random Gaussian weights, or by

Rademacher random variables w_ijε

{−

1, 1}, or by

random phases w_ij=exp(√

{square root over (−

1)}a_ij), a_ijε

[−

π

,π

], or by other randomly drawn weights used in sketching techniques, and based on various types of randomized separable, or, nonseparable source/receiver encodings with random signs, Gaussian weights, random phase rotations, or random time shifts, followed by superposition, or various types of randomized separable, or, nonseparable selections of sources/receivers using uniform random sampling with or without replacement, jitter sampling with or without replacement, or importance sampling with and without replacement.

30. A method according to claim 29, wherein the effects of the subset selection are controlled by either redrawing the subsets or by increasing the cardinality of the subsets after each gradient calculation.

31. A method according to claim 22, wherein a suitable regularization term on the physical model parameters is included as in Equation (11):

32. A method according to claim 31, wherein regularization by transform-domain sparsity promotion via (weighted) one-norm minimization is included.

33. A method according to claim 31, wherein regularization by TV-norm minimization is included.

34. A method according to claim 31, wherein regularization by Sobolev-norm minimization is included.

35. A method according to claim 23, wherein the (weighted) least-squares misfit terms in Equation (4) are replaced by misfit functionals that are robust with respect to outliers in the observed data.

36. A method according to claim 23, wherein the (weighted) least-squares misfit terms in Equation (4) are obtained by measuring the misfit in (weighted) matrix norms.

37. A method according to claim 22, wherein nuisance parameters required to compute updates of the physical model parameters are calculated while performing steps (ii) through (v).

38. A method according to claim 23, wherein the penalty formulation of Equation (4) is solved including uncertainty quantification, comprising generation of random realizations for the physical model parameters from a probability distribution function (14)

π

(m)∝

exp(−

(m−

m*)H<sup>−

1</sup>(m−

m*)), 

 

(14)

where m* is a minimizer of the penalty objective and H is the (Gauss-Newton) Hessian of the penalty objective function.

39. A method according to claim 22 for use in the field of PDE constrained optimization problems selected from the group consisting of process control, medical imaging, non-destructive testing, radar imaging, remote sensing, radiative transfer, and design.

40. A method according to claim 22 for use in a field selected from the group consisting of travel-time tomography, geophysical prospecting, electromagnetic inversion in geophysical prospecting, inversion of gravity data in geophysical prospecting, and wave-equation based inversion in geophysical prospecting, and wherein the physical model is a model of a geological area.

41. A method according to claim 39, wherein the use is in the field of wave-equation based inversion in geophysical prospecting, and wherein the model of a geological area has poro/visco-elastic model parameters.

42. A method according to claim 40, wherein a wave-equation-based migration is conducted.

43. A method according to claim 40, wherein a linearized amplitude-versus-offset/angle/ray-parameter inversion is conducted.

44. A method according to claim 40, wherein a wave-equation based migration velocity analysis is conducted.

45. A method according to claim 40, wherein a full-wave inversion is conducted.