CC BY
81
УДК 550.834, 550.837.21, 550.83.017 Вестник СПбГУ. Сер. 4. Т. 2 (60). 2015. Вып. 3
D. М. Molodtsov1, D. Colombo2, | Yu. V. Roslov], V. N. Troy an1, В. M. Kashtan1
COMPARISON OF STRUCTURAL CONSTRAINTS FOR SEISMIC-MT JOINT INVERSION IN A SUBSALT IMAGING PROBLEM*
1 St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation
2 EXPEC ARC, Saudi Aramco
The goal of this work is to compare three structural constraints, which differ by the amount of prior information they inject in the joint inversion. One of these constraints is the cross-gradient (Gallardo and Meju, 2003) and another two force the gradients to be either parallel or antiparallel, depending on prior information about the sign of velocity-resistivity correlation; one of the constraints also forces the gradients to have identical support. We introduce a flexible framework for 2-D joint inversion, based on Gauss—Newton method, which uses independent meshes for velocity, resistivity and the structural constraint. It is used to compare the considered structural constraints on a synthetic seismic and MT data, simulating a marine subsalt imaging scenario. To study theoretical capabilities of the constraints we first consider a toy "model fusion" example, assuming the true resistivity is known, and then proceed to a more realistic joint inversion experiment. Refs 11. Figs 3.
Keywords: seismic tomography, magnetotellurics, inverse problem.
Д. M. Молодцов1, Д. Коломбо2, Ю. В. Рослое , В. Н. Трояп1, Б. М. Каштан
СРАВНЕНИЕ СТРУКТУРНЫХ ОГРАНИЧЕНИИ ДЛЯ СОВМЕСТНОЙ ИНВЕРСИИ СЕЙСМИЧЕСКИХ И МАГНИТОТЕЛЛУРИЧЕСКИХ ДАННЫХ В ЗАДАЧЕ ИЗУЧЕНИЯ ПОДСОЛЕВОЙ ЧАСТИ РАЗРЕЗА
1 Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
2 EXPEC ARC, Saudi Aramco
Цель данной работы состоит в экспериментальном сравнении трёх структурных ограничений, различающихся характером привносимой в совместную инверсию априорной информации. Первое из этих ограничений — это кросс-градиент (Gallardo and Meju, 2003), два других, в зависимости от априорной информации о знаке корреляции между параметрами, требуют параллельности или антипараллельности градиентов моделей; третье дополнительно требует совпадения носителей градиентов. Предлагается основанная на методе Гаусса—Ньютона гибкая схема двумерной совместной инверсии, использующая независимые сетки для скорости распространения сейсмических волн, удельного электрического сопротивления и структурного ограничения. В численном примере рассматривается инверсия сейсмических и магнитотеллурических синтетических данных для сценария изучения подсолевой части разреза. С целью исследования теоретических возможностей структурных ограничений помимо совместной инверсии рассматривается искусственный пример с точной опорной моделью сопротивления. Библиогр. 11 назв. Ил. 3.
Ключевые слова: сейсмическая томография, магнитотеллурический метод, обратные задачи.
1
Introduction. Joint inversion of magnetotelluric (MT) and seismic traveltime data is useful in sub-basalt, subsalt imaging and other settings with complex geology. One of the key problems in the seismic-electromagnetic joint inversion is to define the constraints linking the velocity and resistivity models. Most often some structural constraints are used,
* This work was supported by St. Petersburg State University research grant 11.38.217.2014 and Saudi Aramco.
e. g., cross-gradients [3] or joint total variation [5]. Since the main goal of the seismic-MT joint inversion is to limit the set of equivalent solutions, one may wish to strengthen the constraints as much as possible using all the available prior information. In this paper we present modifications of the structural constraints, introduced in [7] and [6] and compare them with the well-known cross-gradient constraint and with each other, using a common 2-D Gauss—Newton joint inversion framework. For the numerical experiments we consider a synthetic dataset of marine MT and long-offset seismic, modeling a subsalt imaging scenario at the Red Sea shelf.
Structural operators. Let us define a model of slowness s(x) and a model of resistivity r(x), differentiable in some domain Q. Then one can introduce a structural operator c(Vs, Vr), measuring a local similarity between the structures of two models, such that an integral structural similarity can be enhanced by minimizing its L2-norm:
$c(s,r) = ||c(Vs, Vr)||2 = J[c(Vs(x), Vr(x))]2 dx ^ min . (1)
Q
We will refer to $c(s,r) as a structural functional. Measures of structural similarity exist, that do not have a form of an L2-norm, e. g., Gramian [11]. In some papers, e. g. [3], structural similarity is imposed as a hard constraint via the Lagrange multipliers method, which unfortunately is too computationally complex for large-scale problems. Here we focus on the constraints having a form of the least-squares problem (1) with different choice of the structural operator, which fit well in the same Gauss—Newton inversion framework, and we believe such uniformity is essential for a fair comparison. First consider the cross-gradient operator, introduced by Gallardo and Meju [3]:
c®(Vs, Vr) = |Vs x Vr\,
which is zero if the gradients are collinear or one of them is zero. The structural functional $®(s,r) = ||c®||2 is exactly zero when every level set of the model s(x) corresponds to some level set of the model r(x). For example, this holds for any 1-D models and for any models with a functional relation f (s(x),r(x)) = 0, such that f is differentiable and its coefficients are constant. When it is known a priori if s(x) and r(x) correlate positively or negatively, which is often the case in seismic-MT joint inversion, such information can be used to strengthen the structural constraint. Molodtsov et al. [7] proposed a structural operator
cQ(Vs, Vr, h) = |Vs| |Vr| + hVs ■Vr, (2)
with h = ±1 specified a priori and fixed. An elementary identity holds:
4(Vs, Vr) = cQ(Vs, Vr, +1)cQ(Vs, Vr, -1),
which shows that operator (2) is analogous to the cross-gradient, except that it forces the gradients to be parallel or antiparallel, depending on the value of h. When one of the gradients vanishes, the corresponding Frechet derivative of (2) becomes singular, but we observe that for Gauss—Newton optimization only a minor smoothing of the modulus function is required to handle this problem. It is possible to strengthen the constraint further by requiring the flat areas in two models to coincide; to do so one can consider a structural operator
sr
ce(Vs,Vr,/7.) = -— + /*-—, (3)
|Vs|e |Vr|e
where |x|e = y |x|2 + e2 or |x|e = max{|x|, |e|}. Here we use the first option, which allows simple calculation of the exact Jacobian of (3), while in the work [6] the second one was used. If one of the gradients is zero and the other is not, c®(Vs, Vr,h) = 0, therefore $®(s, r, h) = ||c®||2 =0 only if Vs and Vr are, depending on h, parallel or antiparallel, and their support is identical. The last condition is especially useful for reconstructing piecewise constant models. Parameter e requires rather fine tuning, because for too small e operator (3) becomes highly nonlinear, making Gauss—Newton method ineffective, while for too big e magnitudes of the gradients become significant. The approach we use here, is to start with relatively big value of e and gradually decrease it as Gauss—Newton iterations proceed; it was used by Chan et al. [1] in a similar problem of minimizing smooth approximation of total variation of a 2-D image.
Joint inversion framework for seismic traveltime tomography and magne-totellurics. The 2-D MT forward problem is solved by a modification of the finite-element algorithm of Wannamaker et al. [10], for Frechet derivatives the algorithm, described by Rodi [8], is used. Calculation of seismic traveltimes and rays is based on the work [9]. The inverse problem is solved for hyperbolic transforms [4] of p-wave slowness v-1(x) and logarithmic resistivity ln p(x), to satisfy inequality constraints vp min < vp < vp max and pmin < P < pmax. The models s(x) and r(x) are discretized using independent meshes, giving the vectors s and r (of different dimension); traveltime tomography behaves better on a uniform mesh, while for MT it is reasonable to use a mesh that coarsens with depth and covers a bigger domain to satisfy the boundary conditions. Domain Q, on which the structural operator is defined, is specified in the intersection of the model domains and is discretized with its own mesh. Interpolation from the seismic mesh and the MT mesh to the "structure" mesh is performed by (nonstationary) Gaussian filters, Gs and Gr respectively. The Gaussian kernel width is defined based on the input mesh spacing and expected MT and seismic resolution; by increasing the width we may force a structural operator to act only on the low spatial wavenumbers of the models. All the considered structural functionals are not convex, thus additional regularization is necessary. Accordingly, the joint inversion objective function is:
$(s, r) = w [$T(s) + as$s(s)] + (1 - w) [$Z(r) + ar$r(r)] + a$c(Gss, Grr), (4)
where (s) and (r) are weighted L2 misfits of traveltimes and logarithmic impedances respectively, $s(s) and $r (r) are regularization terms, as > 0 and ar > 0 are regularization parameters, a > 0 is a penalty coefficient, and w G [0,1] is a scaling factor, that controls relative weight of the seismic and MT parts. The regularization terms are L\ (for s) and L2 (for r) norm of the model gradient with anisotropic weighting. We minimize (4) using the damped Gauss—Newton (Levenberg—Marquardt) algorithm with preconditioned CG as the linear solver. To accelerate the convergence we use several heuristics: as iterations proceed, as and ar are decreased; a is either a "big" constant or is increased, as ideally we want $c = 0; w is increased, as MT is sensitive to lower spatial frequencies of the model, than the tomography. These are to a large extent backed by the experiments considered below and are by no means universal.
Synthetic data for the Red Sea model. A 2-D section was extracted from the synthetic 3-D multiparameter model of the Red Sea crust, described in the paper [2]. Fig. 1, a, b show parts of resistivity and velocity models, corresponding to the seismic inversion domain, used to generate the synthetic data; one can see the water layer, the sedimentary section with the layer of evaporites and the basement. TM and TE modes of MT field were
0
g -5
N -10
-15 0 -5
p, fim
v , km/s
fiV-i 6
4 3 2
v , km/s fi i—6
v , km/s
diV-i 6 i
20 30 x, km
Fig. 1. True resistivity (a); true velocity (b); resistivity by standalone MT inversion velocity models: standalone first-arrival tomography (d); model fusion with true resistivity using c® (e); same using €q (f); same using c® (g):
captions show Gauss—Newton iteration number, RMS traveltime error ST, RMS log-impedance error SZ, relative velocity error Ss and relative log-resistivity error Sr
(c);
5
5
calculated for 15 frequencies in 0.001-1 Hz range at 55 ocean-bottom stations with 1 km spacing. MT data were contaminated by white zero-mean Gaussian noise with standard deviation 3% in apparent resistivity and 15 mrad in phase. For seismic an ocean-bottom-cable geometry was modelled, with 400 m spacing between both sources and receivers, each source recorded by all receivers; the maximum offset is 60 km. First-arrival traveltimes were contaminated by 10 ms zero-mean Gaussian noise. Starting models are 10 fim halfspace and constant-gradient velocity. We measure the accuracy of a reconstructed model m by a relative error Sm = ||m — mtrue|l2 / ||mtrue|l2 ' 100%. Standalone first-arrival tomography does not resolve the velocity model below the evaporite layer (Fig. 1, d). To explore theoretical capabilities of the structural operators we perform a model fusion experiment (in terminology of Haber and Holtzman Gazit [5]): minimize (4) with w =1 and r = rtrue. We use prior information that velocity and resistivity are positively correlated to define structural operators (2) and (3). Fig. 1, e, f show that operator (2) performs somewhat better than the cross-gradient, but is unable to capture the bottom of evaporites, since it does not enforce coincident boundaries. Using operator (3) we are able to reconstruct the whole structure even in the shadow zones (Fig. 1, g). Fig. 1, c shows the resistivity model, obtained by standalone MT inversion. Since MT is sensitive to the subsalt conductive layer, it gives a bigger picture compared to the standalone seismic tomography. In the joint inversion we set bigger weight to MT and primarily aim to improve the velocity model. In Fig. 2, a, b one can notice excellent fit between the structures of resistivity and velocity models, which corresponds to successful minimization of the cross-gradient norm (Fig. 3), however instead of
iter. 20, SZ = 4.7%, Sr = 34.6%
20 30 x, km
p, fim 100 10 1
p, fim 100 10 1
p, fim 100 10 1
iter. 20, ST = 19 ms, Si = 15.1% v km/s
6 5
4
3 2
iter. 19, ST = 20 ms, Si = 12.9% vp, km/s6
5
4
3 2
iter. 23, ST = 20 ms, Si = 9.6% vp, kmA
5
4 3 2
x, km
Fig. 2. Joint inversion results. Resistivity models are cropped to the seismic inversion domain. Resistivity (a) and velocity (6), using c%; resistivity (c) and velocity (d), using cq; resistivity (e) and velocity (f), using c®:
captions show Gauss—Newton iteration number, RMS traveltime error ST, RMS log-impedance error SZ, relative velocity error Ss and relative log-resistivity error Sr
109 108 107 106 105 104 103
o
o,
10 20 iteration
0 10 20 iteration
1000 900 800 700 600 500 400 300 200 100 0
o
0 10 20 iteration
Fig. 3. Terms of the objective function (4) as functions of iteration number for model fusion and joint inversion numerical experiments with different choice of the structural operator: adding highly non-quadratic structural functionals to the objective function significantly slows down the convergence; 1 — separate inversion; 2 — r = rtrue; 3 — cq, r = rtrue;
4 — c®, r = rtrue; 5 — c%; 6 — Cq; 7 — c®
the low-velocity layer there is a higher-velocity layer, accordingly velocity and shape of the basement are very inaccurate. Interestingly, cross-gradient performs somewhat better in the joint inversion than in the model fusion, probably because the piecewise constant reference resistivity model zeroes cross-gradient almost everywhere, and the kernel of the Gaussian filter Gr is too narrow to smooth it. Operator c® yields some improvement of both velocity and resistivity models with respect to cq (Fig. 1, c-f), though it is not as apparent as for the model fusion and now we fail to reconstruct the velocity reversal associated with the bottom of evaporites.
Conclusions. We have shown that by using prior information about the sign of velocity-resistivity correlation it is possible, remaining within the framework of the structural approach, to improve significantly the results of seismic-MT joint inversion with the cross-gradient constraint. Taking into account low resolution and typical targets of MT and seismic traveltime tomography, positive velocity-resistivity correlation is often a reasonable assumption. However in case of insufficient prior information a softer constraint such as cross-gradient may be a better option.
Acknowledgments. We thank Saudi Aramco for financial support and permission to publish this work, Saint-Petersburg State University for the research grant 11.38.217.2014 and Ministry of Education and Science of the Russian Federation for financial support of the applied research project RFMEFI57614X0052.
References
1. ChanT.F., ZhouH.M., Chan R. H. Continuation method for total variation denoising problems // SPIE's 1995 International Symposium on Optical Science, Engineering, and Instrumentation, International Society for Optics and Photonics. 1995.
2. Colombo D., McNeice G. W., Sandoval CurielE., Fox A. Full-tensor CSEM and MT for subsalt structural imaging in the Red Sea: Implications for seismic and electromagnetic integration // The Leading Edge. 2013. Vol. 32. P. 436-449.
3. Gallardo L. A., MejuM. A. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data // Geophys. Res. Lett. 2003. Vol. 30, N 13. P. 1658.
4. Habashy T. M., Abubakar A. A general framework for constraint minimization for the inversion of electromagnetic measurements // Progr. Electromagnet. Res. 2004. Vol. 46. P. 265-312.
5. Haber E., Holtzman GazitM. Model fusion and joint inversion // Surv. Geophys. 2013. Vol. 34, N 5. P. 675-695.
6. Molodtsov D. M., TroyanV.N., Roslov Yu. V., ZerilliA. Joint inversion of seismic traveltimes and magnetotelluric data with a directed structural constraint // Geophys. Prosp. 2013. Vol. 61, N 6. P. 1218-1228.
7. Molodtsov D., Kashtan B., Roslov Yu. Joint inversion of seismic and magnetotelluric data with structural constraint based on dot product of image gradient // Exp. Abstr. SEG Annual Meeting. 2011.
8. Rodi W. L. A technique for improving the accuracy of FE solutions for MT data // Geophys. J. Roy. Astr. Soc. 1976. Vol. 44. P. 483-506.
9. Vsemirnova E. A., Roslov Yu. V. Raytracing on irregular mesh // Abstr. V Int. Conf. "Problems of Geocosmos". 2004. P. 68.
10. Wannamaker P. E., StodtJ.A., Rijo L. A. Stable finite element solution for two-dimensional magnetotelluric modeling // Geophys. J. Roy. Astr. Soc. 1987. Vol. 88, N 1. P. 277-296.
11. Zhdanov M. S., Gribenko A., Wilson G. Generalized joint inversion of multimodal geophysical data using Gramian constraints // Geophys. Res. Lett. 2012. Vol. 39, N 9. P. L09301.
References
1. Chan T.F., Zhou H.M., Chan R.H. Continuation method for total variation denoising problems.
SPIE's 1995 International Symposium on Optical Science, Engineering, and Instrumentation, International, Society for Optics and Photonics. 1995.
2. Colombo D., McNeice G.W., Sandoval Curiel E., FoxA. Full-tensor CSEM and MT for subsalt structural imaging in the Red Sea: Implications for seismic and electromagnetic integration. The Leading Edge, 2013, vol. 32, pp.436-449.
3. Gallardo L.A., Meju M.A. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data. Geophys. Res. Lett., 2003, vol. 30, no 13, pp.1658.
4. Habashy T.M., Abubakar A. A general framework for constraint minimization for the inversion of electromagnetic measurements. Progr. Electromagnet. Res., 2004, vol. 46, pp.265-312.
5. Haber E., Holtzman Gazit M. Model fusion and joint inversion. Surv. Geophys., 2013, vol. 34, no 5, pp.675-695.
6. Molodtsov D.M., Troyan V.N., Roslov Yu.V., Zerilli A. Joint inversion of seismic traveltimes and magnetotelluric data with a directed structural constraint. Geophys. Prosp., 2013, vol. 61, no 6, pp.1218-1228.
7. Molodtsov D., Kashtan B., Roslov Yu. Joint inversion of seismic and magnetotelluric data with structural constraint based on dot product of image gradient. Exp. Abstr. SEG Annual Meeting, 2011.
8. Rodi W.L. A technique for improving the accuracy of FE solutions for MT data. Geophys. J. Roy. Astr. Soc., 1976, vol. 44, pp.483-506.
9. Vsemirnova E.A., Roslov Yu.V. Raytracing on irregular mesh. Abstr. V Int. Conf. "Problems of Geocosmos". 2004, pp.68.
10. Wannamaker P.E., Stodt J.A., Rijo L.A. Stable finite element solution for two-dimensional magnetotelluric modeling. Geophys. J. Roy. Astr. Soc., 1987, vol. 88, no 1, pp.277—296.
11. Zhdanov M.S., Gribenko A., Wilson G. Generalized joint inversion of multimodal geophysical data using Gramian constraints. Geophys. Res. Lett., 2012, vol. 39, no 9, pp.L09301.
Статья поступила в редакцию 8 мая 2015 г.
Контактная информация
Molodtsov Dmitry Mikhailovich — post-graduate student; e-mail: molodtsovd@gmail.com
Colombo Daniele — PhD.
Roslov Yuri Viktorovich (1963-2015) — PhD.
Troyan Vladimir Nikolayevich — Dr. Sci., Professor.
Kashtan Boris Markovich — Dr. Sci., Professor.
Молодцов Дмитрий Михайлович — аспирант; e-mail: molodtsovd@gmail.com Коломбо Даниэле — PhD.
Рослов Юрий Викторович (1963—2015) — кандидат физико-математических наук. Троян Владимир Николаевич — доктор физико-математических наук, профессор. Каштан Борис Маркович — доктор физико-математических наук, профессор.
