A constrained-based optimization approach for seismic data recovery - - PowerPoint PPT Presentation

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

A constrained-based optimization approach for seismic data recovery problems

ICASSP 2014 — SS5: Seismic signal processing

M-Q. PHAM, L. DUVAL* IFP Energies nouvelles 1 et 4 av. de Bois-Pr´ eau, 92852 Rueil-Malmaison - France

• C. CHAUX
• Univ. Aix-Marseille, I2M UMR CNRS 7373

39 rue F. Joliot-Curie, 13453 Marseille - France J-C. PESQUET

• Univ. Paris-Est, LIGM UMR CNRS 8049

5 bd Descartes, 77454 Marne-la-Vall´ ee - France

4-9 May 2014

1 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

The fast way

◮ Here: signal = primary; disturbance = multiples ◮ Use of approximate disturbance models (given)

◮ similar to acoustic echo cancellation, speech

dereverberation, pattern matching. . .

◮ imperfect models: adaptation (and regularization) required

◮ Two aims:

◮ to recover highly under-determined signals ◮ to alleviate and fasten hyper-parameter determination

◮ An entwined approach:

◮ pragmatic: with geophysically-sound functions/models ◮ heuristic: explore other suitable penalties/constraints

◮ Focus on the “sparsity” term (signal wavelet coefficients)

2 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Outline

INTRODUCTION CONTEXT Seismic data Problem formulation GENERAL SCHEME MAP estimation Hyperplane constraints RESULTS Synthetic data Real data CONCLUSION

3 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Seismic data propagation & acquisition

• Towed streamer

Hydrophone Primaries + multiple reflection disturbances

4 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Multiple problem formulation

z(n)

• bserved signal

= ¯ s(n)

• multiple

+ ¯ y(n)

• primary

+ b(n)

• noise

where

◮ n: time index (somehow related to underground depth) ◮ z = (z(n))0≤n<N: observed data, combining:

• 1. primary ¯

y = (¯ y(n))0≤n<N (signal of interest, unknown)

• 2. multiples (¯

s(n))0≤n<N (sum of undesired bouncing signals)

◮ b = (b(n))0≤n<N: additive noise realization.

5 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Availability of different approximate models

◮ Several (J) models (or templates) r(n) j

are given

◮ Imperfect in time, amplitude and frequency ◮ Assumption: models linked to ¯

s(n) and corrected throughout time-varying (FIR) adaptive filters ¯ s(n) =

J−1

• j=0

p′+Pj−1

• p=p′

¯ h(n)

j

(p)r(n−p)

j

where

◮ ¯

h(n)

j

: unknown impulse response of the filter corresponding to model j and time n (Pj tap coefficients)

6 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Model examples

100 200 300 400 500 600 700 800 900 1000

¯ s(n) r(n)

1

r(n)

2

Simulated Multiple, 1st model, 2nd model

7 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Model examples

Receiver number Time (s)

1500 1600 1700 1800 1900 1.5 2 2.5 3 3.5 4 4.5 5 5.5

a) Receiver number

1500 1600 1700 1800 1900 1.5 2 2.5 3 3.5 4 4.5 5 5.5

b)

Real data Model

7 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Multiple problem reformulation

z

• bserved signal

= R ¯ h

• filter

+ ¯ y

• primary

+ b

• noise

where

◮ ¯

s = J−1

j=0 Rj¯

hj = R¯ h = ¯ s(0) · · · ¯ s(N−1)⊤

◮ R =

• R0 · · · RJ−1
• and ¯

h =

• ¯

h⊤

0 · · · ¯

h⊤

J−1

⊤

◮ ¯

hj =

• ¯

h(0)

j

(p′) · · · ¯ h(0)

j

(p′ + Pj − 1) · · · · · · ¯ h(N−1)

j

(p′) · · · ¯ h(N−1)

j

(p′ + Pj − 1) ⊤

◮ Rj is a block diagonal matrix.

8 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Maximum-a-Posteriori (MAP) estimation

◮ Assumption: ¯

y is a realization of a random vector Y, whose probability density is given by: (∀y ∈ RN) fY(y) ∝ exp(−ϕ(y))

◮ Assumption: ¯

h is a realization of a random vector H, whose probability density can be expressed as: (∀h ∈ RNP) fH(h) ∝ exp(−ρ(h))

◮ Assumption: b is a realization of a random vector B, of

probability density: (∀b ∈ RN) fB(b) ∝ exp(−ψ(b)) B is assumed to be independent from Y and H

9 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Maximum-a-Posteriori (MAP) estimation

MAP estimation of (y, h) minimize

y∈RN,h∈RNP

ψ

• z − Rh − y
• fidelity: linked to noise

+ ϕ(y)

• a priori on the signal

+ ρ(h)

• a priori on the filters

.

9 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Constraints

Assumption: Gaussian noise, B ∼ N(0, σ2) ψ = ιS where S =

• w ∈ RN | w2 ≤ Nσ2

Assumption: slow variations of adaptive FIR filters along time

• ∀(j, n, p)
• |h(n+1)

j

(p) − h(n)

j

(p)| ≤ εj,p (hyperslab C) Assumption: limits on (functions of) primary coefficients Fy ∈ D, where F ∈ RK×N: analysis frame operator

10 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Constraints

Assumption: Gaussian noise, B ∼ N(0, σ2) ψ = ιS where S =

• w ∈ RN | w2 ≤ Nσ2

Assumption: slow variations of adaptive FIR filters along time

• ∀(j, n, p)
• |h(n+1)

j

(p) − h(n)

j

(p)| ≤ εj,p (hyperslab C) Assumption: limits on (functions of) primary coefficients Fy ∈ D, where F ∈ RK×N: analysis frame operator MAP estimation of (y, h) minimize

y∈RN,h∈RNP ϕ(y) + ρ(h) + ιS

• z − Rh − y
• + ιC(h) + ιD(Fy).

10 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Hard constraints on primaries y: convex set D

F: analysis frame operator, with L subbands (xk = (Fy)k) D = D1 × · · · × DL ∀ ℓ ∈ {1, . . . , L}

• 1. Dℓ = {(xk)k∈Kℓ |

k∈Kℓ φℓ(xk) ≤ βℓ}, where, βℓ ∈ ]0, +∞[,

φℓ : Rcard(Kℓ) → [0, +∞[ is a proper lsc convex function.

• 2. Dℓ = {(xk)k∈Kℓ |

k∈Kℓ φℓ((Lz)k)xk = βℓ}, where

φℓ : R → R, L ∈ RN×N: linear operator, e.g. a Wiener-like filter L = λ1 Diag

• (1 + λ1 + λ2R(0)2)−1 , . . . ,

(1 + λ1 + λ2R(N−1)2)−1 (λ1, λ2) ∈ ]0, +∞[2 , R(n): n-th row of matrix R.

11 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Qualitative results

100 200 300 400 500 600 700 −1.5 −1 −0.5 0.5 1 1.5

Observed signal z (red; σ = 0.01), original y (blue). D is the intersection of two hyperplanes defined from the identity (Id) and the sign functions.

12 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Qualitative results

100 200 300 400 500 600 700 −1.5 −1 −0.5 0.5 1 1.5

Estimated signal ˆ y (magenta), original signal y (blue). D is the intersection of two hyperplanes defined from the identity (Id) and the sign functions.

12 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Qualitative results

100 200 300 400 500 600 700 −0.5 −0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5

Estimated multiples ˆ s (magenta), original multiples s (blue). D is the intersection of two hyperplanes defined from the identity (Id) and the sign functions.

12 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Quantitative results

σ 0.01 0.04 φ α SNRy SNRs α SNRy SNRs 0.4 23.98 15.79 0.9 15.03 9.60 ℓ∞ 0.6 24.48 15.81 0.8 16.24 8.69 ℓ2 0.6 25.59 16.02 0.8 17.84 9.20 ℓ1 0.4 25.98 16.16 0.9 18.19 6.61 sign 0.3 24.43 14.73 0.1 14.75 4.58 Id 0.4 26.19 15.81 0.2 19.74 8.84 Id + sign 0.3 26.40 15.56 0.1 18.43 5.94 Upper table part: “classical constraints” and lower table part: hyperplane contraints. (ϕ = (1 − α)ℓ1, ρ = αℓ2)

13 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Real data

Portion of a common receiver gather: (left) recorded seismic data with a partially appearing primary (right) multiple wavefield template.

14 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Real data

Subtraction results, low field-noise case: primaries (separated from multiples) with (left) [Ventosa et al. 2012] (right) with the proposed method.

14 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Real data

Subtraction results, high field-noise case: primaries (separated from multiples) with (left) [Ventosa et al. 2012] (right) with the proposed method.

14 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Conclusion: in symbols

minimize

y∈RN, h∈RQ

(1 − α)ϕ(y) + αρ(h) subject to      ψ(z − y − Rh) ≤ 1 h ∈ C Fy ∈ D with easier bounds inferred from data & coarser processing

15 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

Conclusion: in words

◮ New variational framework for multiple removal in

seismic data

◮ Use of novel convex proximal splitting algorithms ◮ Generic methodology to impose sparsity and regularity

properties through constrained adaptive filtering in a (wavelet frame) transformed domain

◮ Different strategies permitted for:

◮ sparse modeling, ◮ additive noise removal, ◮ adaptive filter design under appropriate regularity ◮ amplitude concentration constraints

◮ Versatility of the optimization framework

◮ usable in echo cancellation, de-reverberation, matching 16 / 17

INTRODUCTION CONTEXT GENERAL SCHEME RESULTS CONCLUSION

More for free: additional references

Pham, M. Q., C. Chaux, L. Duval, L. and J.-C. Pesquet, A Primal-Dual Proximal Algorithm for Sparse Template-Based Adaptive Filtering: Application to Seismic Multiple Removal: IEEE Transactions Signal Processing, to be published, 2014 http://www-syscom.univ-mlv.fr/˜mpham/publications.html http://arxiv.org/abs/1405.1081 Ventosa, S., S. Le Roy, I. Huard, A. Pica, H. Rabeson, P. Ricarte, and L. Duval, Adaptive multiple subtraction with wavelet-based complex unary Wiener filters: Geophysics, vol. 77, V183–V192, Nov-Dec. 2012 http://arxiv.org/abs/1108.4674 Jacques, L., Duval, L., Chaux, C. and Peyr´ e, G., A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity: Signal Process., vol. 91, 2699–2730, Dec. 2011 http://arxiv.org/abs/1101.5320 Combettes, P. L. and Pesquet, J.-C., Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators: Set-Valued and Variational Analysis,

• vol. 20, no. 2, pp. 307–330, Jun. 2012

http://www.ljll.math.upmc.fr/˜plc/ J.-C. Pesquet and N. Pustelnik, A Parallel Inertial Proximal Optimization Method: Pacific Journal of Optimization, vol. 8, no. 2, pp. 273–305, Apr. 2012 http://perso.ens-lyon.fr/nelly.pustelnik/ 17 / 17