Совместная двумерная инверсия данных электротомографии и РМТ/АМТ

Каминский А.Е.

Zond software kaminae@yandex.ru  

Аннотация

Разработан алгоритм для совместной двумерной интерпретации данных электротомографии и магнитотеллурических зондирований. Алгоритм реализован в программе ZondRes2D и протестирован на ряде синтетических и полевых примеров. Результаты тестирования показали высокую эффективность данного подхода и значительное повышение качества получаемых результатов при решении сложных геологических задач.  

Ключевые слова: электротомография, магнитотеллурические зондирования, совместная инверсия геофизических данных.

Комплексирование геофизических методов один из наиболее эффективных приемов повышения надежности результатов интерпретации [1]. Комплексирование геофизических методов уже достаточно давно существует как отдельное направление разведочной геофизики. Удачное сочетание геофизических методов порой позволяет решить самые сложные геологические задачи.

Вопросами совместной интерпретации данных электрических и электромагнитных зондирований занимаются более 40 лет.  К настоящему времени известно более десятка работ [2,4] в данной области, однако все они не были доведены до программного исполнения. Отсутствие  подобного инструмента интерпретации послужило толчком к разработке программы  для совместной двумерной инверсии данных электротомографии и магнитотеллурических зондирований.

Совместная интерпретация данных электрических и электромагнитных методов является надежным инструментом повышения качества получаемых результатов. Улучшение качества в данном случае, связано с различной чувствительностью электрических и электромагнитных методов к параметрам геоэлектрического разреза. При этом значительно сужается область эквивалентности решения обратной задачи [4].

Электротомография успешно применяется при решении многих геологических задач, однако, обладает рядом недостатков: довольно низкой производительностью, сильным влиянием верхних экранирующих плохо проводящих отложений, а также сравнительно небольшой глубинностью и вертикальным разрешением. Кроме того, существует проблема эквивалентности тонкого слоя, заключающаяся в невозможности раздельного определения его удельного сопротивления и мощности.

В отличие от электротомографии, электромагнитные зондирования характеризуются высокой производительностью и большой глубинностью. Плохо проводящие горизонты не являются препятствием для распространения электромагнитного поля, что позволяет проводить работы в зимнее время и в районах развития высокоомных приповерхностных отложений. С другой стороны, электромагнитные зондирования чаще всего  малопригодны для изучения верхней части разреза на глубинах до нескольких метров и менее чувствительны к слоям высокого сопротивления. Перечисленные выше особенности электрических и электромагнитных методов  наводят на мысль об их взаимной дополнительности и целесообразности их совместного применения при изучении геоэлектрического разреза (Рис.1).

Рис.1. Схема получения результата при совместной инверсии.

В 2011 году нашей группой был представлен  алгоритм для одномерной совместной инверсии данных  электрических и электромагнитных зондирований [3], реализованный в серии одномерных программ Zond.  В данном докладе мы предлагаем двумерное решение данной задачи.  

Нами реализован  алгоритм для совместного решения обратной задачи электротомографии и магнитотеллурических зондирований (в широком частотном диапазоне) на базе программ ZondRes2D [5] и ZondMT2D. Решение использует совместный функционал невязки, контролирующий значимость того или иного типа данных. Для борьбы с гальваническими эффектами в магнитотеллурическом методе использовалось методика подбора статических сдвигов при инверсии.

Главной проблемой при совместной интерпретации данных является разная масштабность комплексируемых методов, то есть разные длины волны по отношению к геологическому объекту и площади интегрирования регистрируемого сигнала. Это может иногда приводить к появлению ложных объектов на результирующем двумерном изображении. 

Рис.2. Результаты тестирования алгоритма на синтетическом примере.

Еще одна проблема возникает при недостаточном глубинном перекрытии электротомографии и магнитотеллурических зондирований. В этом случае в разрезе может появиться своеобразная буферная зона, природа возникновения которой связана с  нехваткой данных того или иного метода.

Рис.3. Результаты совместной 2D инверсии данных электротомографии и аудиомагнитотеллурических зондирований.

Разработанный алгоритм опробован на синтетических (Рис.2) и полевых примерах. На рисунке 3 показан результат совместной интерпретации данных электротомографии и АМТ, любезно предоставленных компанией ”Северо-Запад”.  Каждый из использованных методов дополняет другой, что позволяет лучше разрешать параметры геоэлектрического разреза, т.е. решить геологическую задачу.

Рис.4. Результаты совместной 2D инверсии данных электротомографии и радиомагнитотеллурических зондирований.

На рисунке 4 показан результат совместной интерпретации данных электротомографии и РМТ, любезно предоставленных компанией ”Ленгипротранс”. Использование совместной инверсии позволило существенно улучшить разрешение результатов и подавить эквивалентность второго слоя.

 Полученные результаты свидетельствуют о необходимости совместной интерпретации данных электротомографии и магнитотеллурических зондирований при решении сложных геологических задач.

Литература

1.      Новые подходы к интерпретации геофизических материалов, Каминский А.Е., Лухманов В.Л., Тезисы докладов EAGE, Engineering Geophysics, 2012.

2.      Joint Inversion of Geophysical Data, Vozoff  K., Jupp D.,  Geophysical Journal of the Royal Astronomical Society, Vol. 42, No. 3, 1975

3.      Программа совместной интерпретации данных электрических и электромагнитных зондирований NWCOMPLEX1D, Каминский А.Е., Яковлев Д.В., Алексанова Е.Д. Тезисы докладов EAGE, Engineering Geophysics, 2011.

4.      The joint use of transient electromagnetic and vertical electric sounding in groundwater prospecting, Tarasov V, A Kaminsky A., Extended abstract, Geophysics of the 21st Century - The Leap into the Future, 2003.

International EM induction workshop 2014

International EM induction workshops

Germany 24 – 30 August 2014

IAGA (International Association of Geomagnetism and Aeronomy) is one of eight associations within IUGG (International Union of Geodesy and Geophysics) and IUGG is one of thirty unions within ICSU (International Council for Science).

The main activity of working group I.2 of IAGA is the biennial Electromagnetic Induction in the Earth Workshops, initiated in 1972. The EM induction workshops are the most important events for researchers around the world in the field of electromagnetism in geophysics.

The 22nd EM Induction Workshop in 2014 is organized by the German EM community.

Scientific Program

Please check back in early 2014 for more information on submitting abstracts and the scientific program.

Workshop Invited Reviewers

We are pleased to announce invited reviews from a range of acclaimed international scientists who are experts in their fields. These invited reviews provide a unique opportunity to hear the state-of-the-art in electromagnetic methods and applications. As in previous Workshops, these reviews will later be published.

Joint inversion of refraction and reflection seismics traveltimes data

Joint inversion of refraction and reflection seismics traveltimes data

Kaminsky Alex, Zond software

Introduction

The joint inversion method for  arrival times  of refracted and reflected waves is discussed in this paper. Nowadays this approach is  popular and widely used for solving task of deep seismic exploration. Here we consider an algorithm of forward and inverse problems  specialized for near surface  seismics.

In engineering seismics there was many attempt separation of the reflected waves for extra information. It well known that using of refracted waves only, gives too rough  velocity profile.

Unfortunately, in most cases of engineering studies, there is difficult to identify the reflected waves, because it is "noised" by other types of waves. Since the problem of identification of the reflected waves is an independent, complex direction -  and will not be considered here. Let us assume that for our task, the first breaks and the time of arrival of the reflected waves for the n-th number of boundaries are known.

Forward problem

Forward problem - calculation the first arrival times and the arrival of the reflected waves for a model with an arbitrary velocity distribution and geometry of the reflector based on the Shortest path method (Moser, 1991). This method allows us to calculate the shortest ray path for refracted waves.

The combination of minimum paths  from the source and the receiver to the reflector point allows the build path of the reflected wave for each boundary. Minimal total travel time from the source and the receiver is selected as a reflection point of the border. Using of this method removes the restriction on the complexity geometry  of topography and the  reflecting boundary and  can be used for the interpretation of engineering seismics .

Forward problem has been developed for the two types of models. In the first case velocity section defined as a set of layers with arbitrary geometry boundaries and arbitrary distribution of the velocity inside each layer. The complexity of the boundaries is controlled by nodes number. Any boundary can be reflected and refracted, or only refracted. The advantages of this model type is the possibility of joint interpretation of P and S waves in a common geometry boundaries. Also, it is convenient to use for  sparse observation systems.

In the second case, the model is divided regular mesh of cells with arbitrary velocity. Reflected boundaries are set arbitrarily and are not joined with the geometry of the cells. This type of model is useful for dense observation networks, such as seismic tomography.

Additional extra features of the algorithm includes the possibility of taking account of surface  topography,  anisotropy and the attenuation parameter.

Testing algorithm of forward problem was carried out for a number of analytical solutions and using other well known algorithms (Fig. 1). Shortest path method is based on graph theory, and has controlled accuracy, so when a sufficiently dense subdividing of boundaries could get accuracy less than 0.01 percent.

Figure 1 Calculated reflected (A) and refracted (B) travel times curves and ray paths for complex velocity model.

Inverse problem

For the inverse problem Occam inversion (Constable et all, 1987) procedure was used. The main problem with the joint inversion of velocities and geometry of boundaries is the difference of dimensions for this parameters. This negatively affects the properties of the matrix system. To reduce the dynamic range of the matrix - logarithmic norms of parameters (velocities and  local thicknesses of layers and  apparent velocities) were used.

        Using of thicknesses instead  the depths is allow to avoid the problem of boundaries intersection in the resulting model.  Also, in order to suppress  strong oscillations of geometry,   an additional parameter controls the relative speed of change for velocity and thicknesses was used. Smoothing filter is constructed differently for the two types of models. In the first case, the operator is designed separately for each layer and works in the horizontal direction.  In the second case the common filter is used for smoothing all cell's velocities and additional  - for smoothing  layer's boundary.

Algorithm testing

Inversion algorithm was tested on synthetic data (Fig.2) calculated for several types of velocity models. As a tested model - four-layer cross-section  with arbitrary boundaries and the high-velocity object inside second layer was used.

Node sampling interval for boundary correspond to twice distance between the geophones. Observing system correspond to seismic tomography with shot point at each geophone. Directly before each inversion of synthetic travel times noise component was added.

Testing was performed as follows:

·         In the first stage only first arrival times for the first and second types of models  were inverted. The horizontally layered  model with a constant velocity gradient  was used as start.

·         At  the second stage only reflected times for the first and second types of models  were inverted. The horizontally layered  model with a constant velocity  was used as start.

·         In the third stage, the joint  inversion for refracted and reflected waves data was carried out. In our experience, when working with synthetic data, to achieve acceptable data fitting, there is only  three or four iterations enough.

·         Finally, the inverted models were compared to the original.

Figure 2 Example of refracted and reflected synthetic travel times inversion. A – calculated reflection travel times data for original model C; B – calculated refraction travel times data for original model C; C – four-layers arbitrary velocity model with one reflections boundary(3); D – result of joint   reflected and refracted data inversion; E - result of refracted data inversion; F - result of reflected data inversion. Before every inversion pass, noise was added to the synthetic data.

Figure 2 shows the results of testing the algorithm for arbitrary layered model. This is a typical velocity section at engineering geophysical surveys. The upper part of the model consists of a low-velocity sedimentary rocks. The basis of the section is high-velocity rock formations.

Synthetic seismic survey consisted of 24 geophones spaced at 5 meters and  shots in each geophones point.

The inversion was carried out until the RMS did not reach a specified level of noise. For this example, only 4-5 iterations was enough. The average computation time for each inversion was  about  two minutes.

As seen from the figure 2, the best results were obtained from the joint inversion of reflected and refracted data 2D.  It is practically corresponds to the original model.  With the one hand, using of reflection data allow to get the best rays coverage (for a more reliable determination of the velocity), on the other - it is better to recover  the reflecting boundary geometry.

When the inversion carried out for reflected or refracted data separately, the model recovered considerably rougher. This is especially well seen for inversion of refracted data 2E. The lower boundary  (reflected) shifted by almost 10 meters.

Test results for a number of models have shown significant improvement in the accuracy of recovering parameters in the joint inversion. Algorithm for joint interpretation of refracted and reflected travel times data was  realized in the ZondST2D and is currently  testing with field's data.

Conclusions

Proposed algorithm can be used for processing of shallow seismic data. The joint inversion of reflected and refracted travel times data significantly improves the quality of the resulted velocity sections. Using of refracted and reflected travel times data together  allow to achieve greater depth of investigation.

Reference

1. Moser T.J., Shortest path calculation of seismic rays, Geophysics, vol. 56, no. 1, 1991.

     2. Constable S.,Parker R., Constable C. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. 1987. Geophysics 52, 289.

Software for 2d gravity and magnetic

Разработка алгоритма совместной инверсии времен прихода преломленных и отраженных волн

Каминский А.Е.

http://zond-geo.com

Введение

В работе рассматривается методика совместной инверсии времен прихода преломленных и отраженных волн.  Данное направление сейчас популярно и широко используется  за рубежом для решения задач глубинной сейсморазведки.

Основными  подходами к совместной интерпретации данных преломленных и отраженных волн являются:

• Полноволновая инверсия во временной или частотной области (Zhou et all, 2003)
• Совместная инверсия времен прихода отраженных волн и первых вступлений (Hobro et all., 2003)

Первое направление имеет ряд существенных ограничений при использовании в инженерной сейсморазведке. Второе направление, применительно к малоглубинным изысканиям, рассматривается в рамках данной работы.

В инженерной сейсморазведке уже достаточно давно совершаются попытки выделения отраженных волн для получения дополнительной информации.  Ведь как известно, использование только рефрагированных волн при интерпретации дает слишком грубое представление о скоростном разрезе (рис.1).  

К сожалению, в большинстве случаев, для данных полученных при инженерных изысканиях, выделить отраженные волны затруднительно, так как они сильно “зашумлены” другими типами волн.  Так как задача выделения отраженных волн является самостоятельным,  чрезвычайно сложным направлением -  здесь она рассматриваться не будет.  Примем, что для определенного профиля выделены первые вступления и времена прихода отраженных волн для n-го количества границ.

Рисунок 1 Результаты двумерной инверсии полевых данных сейсмотомографии на рефрагированных волнах.

Предпосылки к совместной интерпретации преломленных и отраженных волн

Основные минусы сейсморазведки на рефрагированных волнах хорошо известны, это:

• Относительно малая глубинность исследований
• Определенные  требования к скоростному разрезу
• Низкая разрешающая способность  

Интерпретация времен прихода тоже имеет ряд ограничений, особенно  при изучении верхней части разреза:

• Отраженные волны часто “загрязнены” другими типами волн
• Трудности корреляции при сложной геометрии границ
• Годографы имеют более сложный вид

Анализ плюсов и минусов обоих направлений приводит к мысли о совместной инверсии времен прихода преломленных и отраженных волн. Это может дать нам следующие преимущества:

• Хорошо разрешение в верхней части дают рефрагированные волны
• Отраженные волны позволяет увеличить глубинность и устойчивость результата за счет хорошего лучевого покрытия (даже при наличии только одной отражающей границы)
• Более четкое определение геометрии отражающей границы за счет правильного определения скоростного разреза

Основным условием совместного использования является не слишком большое различие глубины проникновения преломленных волн и опорной отражающей границы.

Расчет времен первых вступлений и отраженных волн

Прямая задача – расчет времен первых вступлений и  прихода отраженных волн для среды с произвольной распределением скоростей и геометрии отражателей базируется на методе Shortest path (Moser, 1991).  Этот метод позволяет рассчитать кратчайший путь, по которому проходит рефрагированная волна. Комбинация траекторий лучей минимального пробега от источника и приемника к отражателю позволяет построить путь отраженной волны для каждой границы.  В качестве точки отражения выбирается участок границы с минимальным суммарным временем пробега  от источника и приемника. Таким образом алгоритм расчета сводится к трем основным этапам:

• Для каждой пары источник-приемник, рассчитывается два нисходящих фронта к отражающей границе
• Вдоль отражающей границы получаем 2 распределения времен пробега от источника и приемника
• Точка отражения выбирается исходя из двух условий: минимальная сумма времен пробега, угол падения равен углу отражения

Использование данного метода снимает ограничение на сложность рельефа и углы наклона отражающей границы, что позволяет применять его для интерпретации данных инженерной сейсморазведки.

Прямая задача была реализована для двух типов модели.  В первом варианте скоростной разрез задается набором слоев с произвольной геометрией границ и произвольным распределением скорости вдоль профиля в каждом слое (рис. 2). Сложность геометрии границ контролируется количеством узлов. Любая граница может быть отражающей и преломляющей, либо только преломляющей.  К преимуществам данного варианта модели следует отнести возможность совместной интерпретации P и S волн в рамках одной геометрии границ.  Также его удобно использовать при разреженной системе наблюдений. Во втором варианте модель разбита регулярной сетью ячеек с произвольными значениями скоростей.  Отражающие границы задаются произвольно и не связаны с геометрией ячеек. Такой тип модели удобно использовать для плотных сетей наблюдений, таких как сейсмотомография.

Рисунок 2 Годографы и лучевая схема отраженных(A) и преломленных (B) волн для четырехслойной среды с произвольным распределением скоростей и геометрией границ.

К дополнительным положительным характеристикам алгоритма следует также отнести возможность естественного учета рельефа поверхности измерений, анизотропии скоростей и параметра затухания.  

Тестирование алгоритма решения прямой задачи проводилось для ряда аналитических решений и с использованием сторонних алгоритмов.  Shortest path метод основан на теории графов, и обладает контролируемой точностью, поэтому при достаточно плотном разбиении границ удавалось добиться невязки менее 0.01 процента.

Алгоритм решения обратной задачи

Для решения обратной задачи использовался, уже ставший классическим линеаризованный метод наименьших квадратов в модификации Occam (Constable et all, 1987).  Основной  сложностью при совместной инверсии скоростей и положений узлов границ является несоответствие их размерностей. Это крайне негативно сказывается на свойствах конечной системы. Для уменьшения динамического диапазона матрицы, при инверсии были использованы логарифмические нормы параметров – скоростей и локальных мощностей слоев и логарифмические нормы кажущихся скоростей.  Переход от глубины границы слоя к мощности позволил избежать проблемы пересечения границ в результирующей модели. Также, чтобы подавить сильные осцилляции границ, был использован дополнительный параметр контролирующих относительную скорость изменения  скоростей и границ.  Сглаживающий оператор строился по-разному для двух типов моделей. В первом случае сглаживание производится отдельно для каждого слоя в горизонтальном направлении. Во втором варианте строится общий оператор для сглаживания скоростей в вертикальном и горизонтальном направлении и дополнительный для сглаживания мощностей слоев для соседних узлов границы внутри слоя.

Тестирование инверсии

Алгоритм инверсии тестировался на синтетических данных, полученных для нескольких типов скоростных разрезов.  В качестве основного, был выбран четырехслойный разрез с произвольной геометрией границ, наличием аномальных по скорости объектов внутри слоев и отражающем основанием.  Для основания разреза скорость упругих волн была постоянна. Шаг дискретизации границы соответствовал удвоенному расстоянию между сейсмоприемниками. Система наблюдений соответствовала сейсмической томографии с пунктом возбуждения на каждом сейсмоприемнике. Непосредственно, перед каждым циклом инверсии на синтетические  годографы накладывалась шумовая составляющая разной амплитуды.  

Рисунок 3 Пример инверсии синтетических данных отраженных и преломленных волн. A – расчетные годографы отраженных волн для модели C; B – расчетные годографы преломленных волн для модели C; C – оригинальная модель с отражающим основанием; D – результат совместной инверсии отраженных и преломленных волн; E - результат инверсии только преломленных волн; F - результат инверсии только отраженных волн.

   

Рисунок 2 демонстрирует результаты тестирования алгоритма для произвольной слоистой модели. Верхняя часть разреза представлена низкоскоростными осадочными образованиями, в основании лежат плотные скальные породы.

Система наблюдений состояла из 24 геофонов расположенных через 5 метров, с пунктом возбуждения на каждом.

Как видно из рисунка лучший результат получен при совместной инверсии данных отраженных и преломленных волн.

Тестирование производилось следующим образом:

На первом этапе инвертировались только времена первых вступлений для первого и второго типов модели. В качестве начального приближения использовалась горизонтально-слоистая среда с постоянным градиентом скоростей.

 На втором этапе скоростной разрез подбирался  по временам прихода отраженных волн для первого и второго типов модели. В качестве начального приближения использовалась горизонтально-слоистая среда с постоянным значением скорости.

На третьем  этапе проводилась совместная инверсия данных преломленных и отраженных волн. Как показал наш опыт, при работе с синтетическими данными, для достижения приемлемой невязки достаточно было трех-четырех итераций.

В заключение, полученные модели сравнивались с исходной.

Результаты тестирования для широкого класса моделей показали существенное улучшение точности определения параметров при совместной инверсии данных.

Также было проведено исследование влияния количества данных (рис.4) и шумовой составляющей (рис.5) на результат инверсии.

Рисунок 4  Пример инверсии синтетических данных отраженных и преломленных волн при полной и разреженных системах. A – 24 пункта возбуждения; B – 12 пунктов возбуждения; C – 6 пунктов возбуждения.

Анализ  результатов для систем наблюдений с различной плотностью пунктов возбуждения показал не слишком сильное снижение разрешающей способности.

Рисунок 5  Пример инверсии синтетических данных отраженных и преломленных волн при разном уровне шумовой составляющей.  A – оригинальная модель;. B – восстановленная, шум – 5 мсек; C - шум –10 мсек.  

Алгоритм  совместной интерпретации  данных преломленных и отраженных волн реализован в программе ZondST2D  и в настоящее время проходит стадию тестирования на полевых материалах (рис.6).

Рисунок 6  Результат инверсии полевых данных (отражающая граница в основании разреза).  A – отраженные и преломленные;. B –только преломленные.

Для тестирования на полевых материалах использовалась система наблюдений с 48 геофонами через 2 метра и 9 пунктами возбуждения через 12 метров. Количество данных участвующих  в инверсии составляло: 362 - для рефрагированных,  160 - для отраженных волн. Конечная относительная невязка после инверсии  составила 5.2 процента.

 

Выводы

Предложенный в данной работе алгоритм может быть использован для обработки данных, как “большой”, так и  малоглубинной сейсморазведки. Совместная инверсия данных отраженных и преломленных волн существенно повышает качество результирующих скоростных разрезов.

Ссылки

1. Zhou B., Greenhalgh, S.A., Crosshole seismic inversion with normalized full-waveform amplitude data, Geophysics, 68, 2003.
2. J.W.D.Hobro, S.C.Singh,T.A.Minshull, Three-dimensional tomographic inversion of combined reflection and refraction seismic traveltime data. Geophysical Journal International, 152, 2003.
3. Moser T.J., Shortest path calculation of seismic rays, Geophysics, 56, 1991.
4. Constable S.,Parker R., Constable C. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics 52, 1987.