CC BY
7_ДОКЛАДЫ АН ВШ РФ_
2018_октябрь-декабрь_№ 4 (41)
- ТЕХНИЧЕСКИЕ НАУКИ -
УДК 550.837:517.958
О ПАРАМЕТРИЗАЦИИ ГЕОЭЛЕКТРИЧЕСКОЙ МОДЕЛИ В ЗАДАЧАХ АЭРОЭЛЕКТРОРАЗВЕДКИ В СРЕДАХ С РЕЛЬЕФОМ И СЛОЯМИ ПЕРЕМЕННОЙ ТОЛЩИНЫ
Д.С. Киселев, Н.В. Кондратьев, Ю.И. Кошкина, Д.В. Вагин, М.Г. Персова, Ю.Г. Соловейчик
Новосибирский государственный технический университет
Работа посвящена методу геометрической 3Б-инверсии данных аэроэлектроразведки во временной области в существенно неоднородных геологических средах, характеризующихся резкими перепадами высот рельефа поверхности Земли, изменяющимися толщинами контрастных по проводимости слоев, наличием латеральных неоднородностей проводимости в слоях, перекрывающих локальные целевые объекты. В рассматриваемых ситуациях, во-первых, невозможно «визуально» по форме снятых сигналов определить наличие и положение целевого объекта, а во-вторых, неучет сложной формы рельефа поверхности Земли и границ между другими слоями геоэлектрической модели может приводить либо к пропуску целевых объектов, либо к появлению ложных аномалий. Метод основан на использовании специальной параметризации, которая включает в себя физические свойства структурных частей геоэлектрической модели и их геометрические характеристики. Важной отличительной особенностью является параметризация изогнутых границ между слоями геоэлектрической модели с помощью опорных точек, определяющих плавную деформацию этих границ в процессе нелинейной 3Б-инверсии. Для описания таких поверхностей используются бикубические сплайны. Работоспособность предлагаемого подхода показана на примере обработки синтетических данных, полученных для сложной геоэлектрической модели среды с рельефом и верхним проводящим латерально неоднородным слоем переменной толщины, построенной как обобщение результатов интерпретации данных электромагнитных зондирований по ряду площадей Восточной Сибири. Показано, что за небольшое количество итераций местоположение (как в плане, так и по глубине) целевого объекта в виде кимберлитовой трубки и его основные характеристики были определены с достаточно высокой точностью.
Ключевые слова: аэроэлектроразведка, метод конечных элементов, 3Б-инверсия, рельеф, электромагнитное поле.
Б01: 10.17212/1727-2769-2018-4-77-92
Введение
При проведении аэроэлектромагнитных исследований важнейшую роль играет качество обработки получаемых данных. Одной из основных задач аэроэлектроразведки является поиск твердых полезных ископаемых на небольших и средних глубинах [1, 2]. Целевые объекты (рудные тела, кимберлитовые трубки) в этом случае имеют небольшие размеры, а вмещающие их геологические среды характеризуются, как правило, латерально неоднородными по проводимости и по тол-
Работа выполнена при финансовой поддержке гранта Президента Российской Федерации для государственной поддержки молодых российских ученых - кандидатов наук (№ гранта МК-5010.2018.9).
© 2018 Д.С. Киселев, Н.В. Кондратьев, Ю.И. Кошкина, Д.В. Вагин, М.Г. Персова, Ю. Г. Соловейчик
щине верхними слоями, перекрывающими слои с целевыми объектами, и довольно существенными перепадами высот в рельефе поверхности Земли.
Все эти факторы делают практически бессмысленным применение распространенных на сегодняшний день локальных одномерных инверсий для определения местоположения целевых объектов. В настоящее время опубликовано достаточно много работ, посвященных вопросам 3Б-инверсии. Многие авторы рассматривают мелкоячеистые инверсии (основанные на разбиении исследуемого объема на мелкие ячейки), в которых определяются значения удельной электрической проводимости [3-7]. Однако эти инверсии, в зависимости от уровня регу-ляризирующих добавок в минимизируемом функционале невязки между расчетными и практическими данными, дают либо слишком «гладкие» геоэлектрические модели, либо, наоборот, лишком «пестрые» картины проводимости. В первом случае как минимум очень плохо определяются границы целевых объектов, а как максимум целевые объекты могут быть вообще не выявлены на фоне других «размытых» неоднородностей. Во втором же случае резко растет количество «эквивалентных» решений и вполне могут быть получены ложные целевые объекты. По сути, эти манипуляции, нацеленные на фактически «ручную» настройку параметров регуляризации, делают систему инверсии субъективной, зависящей от конкретного пользователя, что резко снижает доверие к получаемым геологическим моделям.
Альтернативным подходом является 3Б-инверсия, основанная на представлении геоэлектрической модели геологической среды с использованием комбинации физических (проводимости) и геометрических (описывающих границы структурных частей геоэлектрической модели) параметров [8-14]. Применение таких подходов позволяет существенно снизить число искомых параметров геоэлектрической модели, тем самым снижая вычислительные затраты и число возможных эквивалентных моделей и соответственно улучшая сходимость инверсии. Кроме того, это позволяет практически без увеличения вычислительных затрат использовать адаптивные регуляризации, что освобождает пользователя от необходимости каких-либо действий по управлению процессом нелинейной инверсии, улучшает сходимость инверсии и обеспечивает возможность получения «геологически» адекватных моделей.
Однако в работах [9-13], в которых были использованы геометрические 3Б-инверсии, не рассматриваются такие часто встречающиеся на практике ситуации, когда рельеф поверхности Земли и искривленные границы между геологическими слоями характеризуются резкими перепадами высот, в то время как в работах [15, 16] показано существенное влияние этого фактора на измеряемый сигнал. В работе [17] был рассмотрен подход к геометрической 3Б-инверсии для ситуации с искривленными, но фиксированными геологическими границами, что несколько сужает ее возможности.
В данной работе мы рассмотрим метод 3Б-инверсии, также основанный на совместном поиске геометрических и физических параметров, но в котором в искомые параметры геоэлектрической модели будут входить не только параметры, описывающие границы восстанавливаемых неоднородностей по латерали внутри геологических слоев и значения проводимости в них, но и параметры, характеризующие форму искривленных границ между слоями, а также проводимости слоев вмещающей среды. Такая параметризация позволит наиболее полно описать геологическую модель, снизить зависимость от априорных данных, и в результате повысить достоверность восстанавливаемых геоэлектрических моделей геологической среды.
1. Параметризация геоэлектрической модели геологической среды
в задачах аэроэлектроразведки и математический аппарат ЭБ-инверсии
Геоэлектрическая модель представляется набором субгоризонтальных слоев, количество которых может быть взято из априорных геологических данных или определено из Ш-инверсии. В каждом из слоев в рамках обследуемой области задается блочная структура, представляющая собой ряды блоков, ориентированных вдоль некоторого направления. Будем считать, что это направление совпадает с осью X. Тогда неизвестными параметрами могут являться х -координаты границ отдельных блоков, у -координаты рядов блоков, а также значения удельной электрической проводимости внутри каждого блока. Каждая из внутренних границ как между соседними блоками, так и между рядами блоков характеризуется одним параметром, что позволяет восстанавливать геоэлектрическую модель как «непрерывную». Внешние (в плане) границы блоков как по х -, так и по у -координате могут перемещаться независимо от внешних границ соседних блоков. Очевидно, что в этом случае при сжатии блоковой структуры (т. е. при смещении ее крайних левых границ вправо или крайних правых границ влево) освободившееся место занимает «внешняя» среда - в данном случае среда с удельной проводимостью соответствующего слоя, значение которой также входит в вектор искомых параметров. Форма верхних и нижних границ блоков (по вертикали, т. е. в определенном смысле «по 2») повторяет форму поверхностей, определяющих границы между соответствующими слоями геоэлектрической модели.
Все поверхности (поверхность Земли и границы между слоями) описываются с помощью бикубических сплайнов [18], которые, в свою очередь, строятся по
наборам точек |(х™, у', г'), р = 1—Рт |, где т - номер соответствующей поверхности (поверхность Земли, границы между слоями); Рт - число точек, использованных для ее описания [19].
При этом предусмотрено два способа параметризации этих поверхностей при решении обратной задачи. В первом случае форма поверхности является фиксированной: например, повторяет рельеф поверхности Земли, которая фиксируется в процессе аэросъемки и включается в исходные данные при выполнении инверсии, или является строго горизонтальной и др. В этом случае параметром является
приращение вдоль вертикальной оси (2) всех точек (х', у"т, ггр ), описывающих
эту поверхность (по которым строится сплайн для представления ее в конечно-элементной сетке [19]). Во втором случае форма поверхности определяется в процессе инверсии. В этом случае задается начальная форма этой поверхности (например, повторяющая рельеф поверхности Земли или строго горизонтальная) и вводится несколько опорных точек, расположенных в плане по некоторой регулярной сетке, с помощью которых эта форма будет изменяться. Параметрами в этом случае будут являться приращения 2-координаты этих точек. Коррекция поверхности осуществляется следующим образом. После нахождения значений этих параметров на очередной итерации нелинейной инверсии с помощью билинейного интерполянта осуществляется расчет приращений значений 2-координат для
всех точек (х', у', г'), по которым строится описывающий подбираемую поверхность сплайн, и строится новая поверхность.
В зависимости от выбранного способа параметризации поверхности (или поверхностей) в вектор искомых параметров для каждой восстанавливаемой границы
геоэлектрической модели добавляется либо один параметр, либо несколько - по числу заданных опорных точек для подбора формы поверхности.
Компоненты вектора параметров, описывающих геоэлектрическую модель, ищутся путем минимизации следующего функционала:
к ь м _
фа (ь) = ЕЕ(с°г/деп (Ь)) (ъш - К) ^ щи, (1)
1=11=1 т=1 ь
где К - количество положений приемно-генераторной установки аэросистемы; Ь - количество временных каналов; 5ег/ (Ь) - разность теоретических ец и измеренных ей значений сигналов ЭДС в приемнике для 1 -го положения приемно-генераторной установки в I -й момент времени; М - количество искомых параметров, являющихся компонентами вектора Ь; Ът - значение т -го параметра Ът на предыдущей итерации.
Значения весовых функций выбираются равными 1/ ей , при этом, если ей принимает значение ниже некоторого заданного порога щ, определяемого уровнем ошибки измерений для каждой конкретной аэросистемы, то значение весовой функции берется равным 1/ е^ . Регуляризирующие параметры ат выбираются адаптивно в процессе решения обратной задачи. Регуляризирующая добавка с параметрами ат обеспечивает поиск очередного приближения искомых параметров в заданных диапазонах значений, что позволяет не нарушать геометрию структурных частей геологической модели и не приводит к нефизическим значениям их удельной электрической проводимости.
Представляя значения дец (Ь) в виде первых двух слагаемых ряда Тейлора в
окрестности значений параметров на предыдущей итерации Ь, получим
8е, (Ь) *8е, (Ь) д(де'1)
_АЪт, (2)
ь=ь
ЛЛ т=1 иит
где АЪт = Ът - Ът .
С б д(де,1)
Способ вычисления производных -— существенно зависит от типа пара-
дЪт
метров.
Если параметр Ът соответствует координате границы между блоками (рядами блоков) или значению проводимости внутри блока, для вычисления производных
д(де'г) решается следующее векторное дифференциальное уравнение:
дЪт
1 - I дА ^ -
— УхУхА+ стт ^^ = (стт -ст)Е . (3)
цо д
В уравнении (3) поле Е является результатом решения прямой задачи, описанной в работе [19], для значений параметров Ът, полученных на предыдущей итерации и определяющих текущее распределение удельной проводимости ст(х, у, 2). Значение функции стт (х, у, 2) отлично от значений функции
ст(х, у, 2) только в местах, где изменилась проводимость из-за приращения значения Ьт параметра Ьт на величину дЬт .
В этом случае значения производных невязок сигналов (ЭДС) в приемниках
д(8р..,) др.., - „т
определяются по значениям А 1 (х, у, 2, /) в виде
дЬт дЬт
д(УхА„т) кг д/ дЬт
где кг - момент приемника.
Вариационная постановка и конечноэлементная аппроксимация для уравнения (3) строится так же, как и при решении прямой задачи [19]. Принцип построения конечноэлементной сетки также является аналогичным [19, 20], однако размеры шагов сетки берутся в два раза большими и увеличивается коэффициент разрядки. В результате вычислительные затраты на решение уравнения (3) по сравнению с решением прямой задачи сокращаются более чем на порядок.
Если параметр Ьт соответствует геометрическим параметрам, описывающим поверхности между слоями геоэлектрической модели, или значениям проводимо-
д(5в,7)
сти слоев, для вычисления производных-— решается задача, полностью ана-
дЬт
логичная прямой задаче [19], но при этом распределение проводимости определяется с учетом изменения соответствующего параметра Ьт на дЬт . В этом случае
значения производных невязок сигналов (ЭДС) в приемниках —— определяются
дЬт
в виде Ух(Ёт -Ё)——, где Ёт является результатом решения прямой задачи в
дЬт
области, полученной в результате изменения значения Ьт параметра Ьт на значение Ьт + дЬт. Заметим, что конечноэлементная сетка для вычисления Ёт должна быть полностью топологически эквивалентна сетке, используемой для решения прямой задачи, т.е. прообразы этих сеток (в шаблонных координатах [19]) должны быть одинаковыми.
Подставляя соотношение (2) в (1), получаем:
К ь
( (
Фа(Ь) = ££ Щ1 Ьп (ь)+£
м
д(8ег/)
V т=1 дЬт
_ДЬт
ь=ь ))
м 2
+ £ат (ДЬт) ^ т^П. (4)
т=1 АЬ
Минимизация (4) приводит к СЛАУ вида
(А + а)ДЬ = f , (5)
где компоненты матрицы А и вектора правой части f вычисляются с помощью соотношений
й й 11 —ьр дь„
/р = -££8ей(Ь)^, Р,* = 1..М . (6)
г=11=1 дЪР
В СЛАУ (5) матрица а является диагональной матрицей с элементами ат на главной диагонали.
Выбор коэффициентов ат на каждой итерации нелинейной инверсии осуществляется следующим образом [11]. Сначала берутся некоторые (достаточно малые) начальные значения ат и определяются значения приращений ДЪт из решения системы (5). Далее проверяется, чтобы каждое из значений Ът = Ът + ДЪт не выходило за границы, установленные для данного параметра. Кроме того, для геометрических параметров, описывающих латеральные границы блоков, контролируется, что значения координат левых границ блоков (и рядов блоков) не превышают значений координат их правых границ и чрезмерно не приближаются к ним. Для геометрических параметров, описывающих поверхности между слоями, контролируется, что 2-координаты точек, по которым строится сплайн соответствующей поверхности, не превышают 2-координаты точек соседней (сверху или снизу) поверхности и чрезмерно не приближаются к ним. Для тех параметров Ът , для которых контролируемое условие нарушается, соответствующие параметры регуляризации ат увеличиваются в к раз (мы использовали к = 2 ) и решение системы (5) повторяется. Поскольку размер системы (5) достаточно мал (несколько сотен), эта процедура практически не требует вычислительных затрат, но позволяет получить новые параметры Ът как минимизирующие
линеаризованный функционал Ф(Ь), так и удовлетворяющие наложенным на Ът ограничениям. В результате сходимость процедуры нелинейной инверсии значительно улучшается.
Обратим внимание на один важный аспект. Для выполнения геометрической инверсии мы формируем равномерную сетку (далее мы будем называть ее «структурной» сеткой), покрывающую область, в которой должна восстанавливаться геоэлектрическая структура геологической среды. Шаг в этой сетке определяется точностью, требуемой для определения координат границ аномальных тел (например, при решении обратных задач аэроэлектроразведки мы берем этот шаг равным 10 м). Блоки стартовой модели конструируются так, чтобы координаты их границ совпадали с координатами этой «структурной» сетки. Приращения ДЪт, найденные из системы (5) на очередном шаге нелинейной инверсии и соответствующие геометрическим параметрам, корректируются таким образом, чтобы очередное значение параметра Ът (координата границы блока или ряда блоков) совпало с ближайшей координатой «структурной» сетки. Заметим также, что при-
д(5е,,) _
ращения дЪт геометрических параметров для расчета производных - бе-
дЪт
рутся равными одному шагу «структурной» сетки.
Использование «структурной» сетки позволяет исключить появление очень мелких и «лишних» шагов в конечноэлементной сетке, что в свою очередь позволяет избежать резкого ухудшения обусловленности матрицы конечноэлементной системы и роста ее размера.
Однако в результате такой коррекции приращений геометрических параметров, найденных из системы (5), нарушается соответствие их значений найденным
значениям приращений физических параметров (т. е. соответствующих проводи-мостей блоков). Поэтому значения Ьт, соответствующие проводимостям, пере-считываются из минимизации функционала (4) при фиксированных (новых) значениях параметров Ь^ = Ь^ + ДЬ^, соответствующих геометрическим параметрам
(т. е. фактически решается подсистема системы (5) с исключенными из нее неизвестными, соответствующими уже найденным и зафиксированным значениям геометрических параметров).
2. Результаты численных экспериментов
Численные эксперименты проводились на синтетических данных, получаемых для геоэлектрической модели, типичной для условий Восточной Сибири. Эта модель была построена как обобщенная по результатам интерпретации аэроэлектро-разведочных работ. Она содержит четыре слоя переменной толщины. В первом (верхнем) слое удельное сопротивление изменяется в диапазоне от 20 Ом • м до 50 Ом-м. Во второй слабопроводящий слой с сопротивлением р = 200 Ом • м помещен локальный объект с размерами 100 х 100 м2 в плане и пониженным до 10 Ом • м сопротивлением. Он имитирует целевой объект - кимберлитовую трубку. Определение этого объекта по данным электромагнитной съемки затруднено следующими факторами.
Верхний слой (перекрывающий целевой) со средним значением р = 45 Ом • м содержит значительные латеральные неоднородности пониженного (до 20...35 Ом • м) сопротивления, а изменения его толщины достигают 60 м (от 64 до 3 м). Помимо этого рассматриваемая геоэлектрическая модель характеризуется большими перепадами высот рельефа поверхности Земли, которые достигают порядка 100 м при смещении на 1 км по латерали. Третий слой геоэлектрической модели с р = 30 Ом • м имеет слабоменяющуюся толщину и расположен на глубинах приблизительно от 240 до 340 м. Ниже этого слоя задана среда с р = 100 Ом • м. Трехмерный вид этой геоэлектрической модели, рассеченной плоскостями у = -4450 м, у = -4250 м, у = -4150 м (это сечение проходит через целевой объект), у = -4050 м, показан на рис. 1.
С использованием 3Б-моделирования для рассмотренной геоэлектрической модели были синтезированы данные аэроэлектромагнитной съемки для 26 временных каналов в диапазоне от 10 мкс до 1 мс вдоль десяти полетных линий (профилей), удаленных друг от друга примерно на 50 м. Общее количество точек измерений составило около 1300. На рис. 2 показаны электрограммы для пяти временных каналов вдоль профиля у и -4150 м (который проходит непосредственно над целевым объектом - кимберлитовой трубкой) и профиля у и -4050 м (расположенного сбоку от трубки).
Из представленных результатов видно, что по виду сигналов невозможно определить даже само наличие трубки (не говоря уже о ее местоположении) при том, что ее влияние в сигнале достаточно существенно и достигает 30 % (заметим, что в рассматриваемых условиях такой уровень сигнала является уверенно измеримым с использованием современной измерительной аппаратуры). Поэтому определить наличие и местоположение целевого объекта в таких условиях становится возможным только с использованием аппарата высокоразрешающей 3Б-инверсии.
в г
Рис. 1 - Трехмерный вид истинной геоэлектрической модели, рассеченной
вдоль плоскостей:
y = -4450 м (а), y = -4250 м (б), y = -4150 м (в), y = -4050 м (г)
Fig. 1 - A 3-D view of the true geoelectrical model cut along the planes:
y = -4450 m (a), y = -4250 m (b), y = -4150 m (c), y = -4050 m (d)
„t=0.08 4=0.12 7 мс 4 мс
.........' 4=0..19 4=0.32 5 мс 8 мс
,..t=0,58 3 мс
........... .......
_ ..... .....t=0. •... t=0. 08.7. мс .24. м.с..
• У' .....t=0. .....t=0. 1.9.5..м.с. '2.8. мс.
........ .....t=0. 583 мс
3200. 3400. 3600. 3800. X, м а
3200. 3400. 3600. 3800. X, м б
Рис. 2 - Электрограммы для пяти временных каналов вдоль профилей y = -4150 м (а) (над трубкой) и y = -4050 м (б) (сбоку от трубки), рассчитанные для истинной модели («практические» данные для инверсии)
Fig. 2 - Data in five time channels along the profiles y = -4150 m (a) (which runs directly over the tube) and y = -4050 m (b) (which runs outside the tube) calculated for the true model (synthetic data for inversion)
Для параметризации геоэлектрической модели были заданы две блочные структуры: в верхнем слое блочная структура содержала 7 х 5 блоков, а блочная структура во втором слое содержала 6 х 4 блоков. Искомыми параметрами были
взяты значения удельной электрической проводимости внутри блоков, х-коор-динаты границ блоков и у-координаты рядов блоков. Эти блочные структуры показаны на рис. 3.
в г
Рис. 3 - Параметризация геоэлектрической модели:
а и б - блочная структура в верхнем слое; в и г - блочная структура во втором слое; а и в - вид сверху; б и г - трехмерный вид
Fig. 3. Parameterization of the geoelectrical model:
(a, b) - a block structure in the top layer; (c, d) - a block structure in the second layer; (a, c) - top views; (b, d) 3-D views
Для поиска формы границы между верхним и следующим за ним (вторым) слоем была использована следующая параметризация. По сети 250 х 200 м (в пределах от 2750 до 4250 м по оси x и от -4500 до -3700 по у) было расставлено 35 опорных точек, приращения которых вдоль вертикальной оси являются параметрами, определяющими изменение формы поверхности между этими слоями. При нулевых значениях этих параметров форма этой параметризуемой поверхности повторяет форму рельефа поверхности Земли на глубине 50 м.
Верхняя и нижняя границы третьего слоя искались строго горизонтальными и в качестве определяющих их параметров были взяты значения z-координат соответствующих плоскостей. Начальные (для инверсии) значения этих z-координат были взяты равными -250 м и -300 м (в то время как средние значения z-координат верхней и нижней границы третьего слоя в истинной модели составляют 240 и 340 м соответственно).
Также в качестве параметров геоэлектрической модели были взяты удельные электрические сопротивления слоев.
Таким образом, общее количество параметров, которые должны быть определены в ходе 3D-инверсии, было 200.
В качестве стартовой геоэлектрической модели была взята однородная среда (с «правильным» рельефом) с удельным электрическим сопротивлением, равным 100 Ом • м.
В ходе нелинейной 3Б-инверсии было сделано 9 итераций, при этом значение
к ь 2
5вй (Ь))
функционала невязки Ф(Ь) = \ 1 =1- изменилось от 0,41 (для стар-
1 КЬ
товой модели) до 0,27.
На рис. 4 представлен трехмерный вид геоэлектрической модели, полученной в результате девяти итераций нелинейной 3Б-инверсии. Эта модель, так же как и истинная модель (см. рис. 1), представлена с рассечением плоскостями у = -4450 м, у = -4250 м, у = -4150 м (через целевой объект), у = -4050 м.
в г
Рис. 4 - Трехмерный вид полученной после 9-й итерации геоэлектрической модели, рассеченной вдоль плоскостей:
y = -4450 м (а), y = -4250 м (б), y = -4150 м (в), y = -4050 м (г) Fig. 4 - A 3-D view of the geoelectrical model obtained after the 9th iteration; this model is cut along the planes: y = -4450 m (a), y = -4250 m (b), y = -4150 m (c), y = -4050 m (d)
Сравнивая полученную в результате 3D-инверсии геоэлектрическую модель с истинной моделью, можно сделать следующие выводы. Несмотря на очень существенно искажающие сигнал факторы, связанные с резкими перепадами высот в рельефе поверхности Земли, с неровным по толщине верхним слоем и неоднородной по латерали его проводимостью, положение целевого объекта и контраст его удельного сопротивления по отношению к удельному сопротивлению содержащего его слоя были определены с высокой степенью достоверности (рис. 1, в и 4, в). Хотя целевой объект в восстановленной модели имеет чуть большие размеры по латерали и чуть более высокое значение удельного сопротивления по сравнению с объектом истинной модели, точка бурения будет определена достоверно и попадет практически в центр искомого объекта.
Заметим, что с достаточно высокой точностью определено не только местоположение в плане целевого объекта, но и его положение по глубине: определенная в результате 3D-инверсии форма поверхности между верхним и следующим за ним (целевым) слоем достаточно хорошо соответствует форме соответствующей поверхности в истинной модели.
Заметим также, что положение верхней и нижней кромки, а также сопротивление третьего слоя также были определены достаточно хорошо.
На рис. 5 представлены электрограммы, рассчитанные для геоэлектрических моделей, полученных после первой итерации 3Б-инверсии (эти кривые показаны светлым цветом) и после последней, девятой итерации (более темные линии). Для сравнения на этом рисунке точками показаны «практические» данные (синтезированные для истинной модели). Из представленных результатов видно, что полученные кривые, рассчитанные для восстановленной в результате инверсии геоэлектрической модели, достаточно хорошо соответствуют «практическим» данным.
________t=o --------1=0 .087мс. .124 мс.
—t=o —t=o ..1.9.5. .мс. .328мс
________t=0 .583 мс
3200.
3400.
3600. 3800. X, м
3200. 3400. 3600.
3800.0. X, м
а б
Рис. 5 - Электрограммы для пяти временных каналов вдоль профилей y = -4150 м (а) (над трубкой) и y = -4050 м (б) (сбоку от трубки), рассчитанные для истинной модели (черные точки), и для геоэлектрических моделей, полученных после первой (светлая кривая) и девятой (темная кривая) итераций 3Б-инверсии
Fig. 5. Data in five time channels along the profiles y = -4150 m (a) (which runs directly over the tube) and y = -4050 m (b) (which runs outside the tube) calculated for the true model (black points) and for the geoelectrical models obtained after the first (light curves) and ninth (dark curves) iterations of the 3-D inversion
Заключение
Представленный в работе программно-математический аппарат нелинейной геометрической 3Б-инверсии аэроданных и предложенный подход к параметризации геоэлектрической модели позволяют с достаточно высокой точностью определять местоположение (как в плане, так и по глубине) целевого объекта в виде кимберлитовой трубки и его основные характеристики в сложных геоэлектрических условиях: резких перепадов высот в рельефе поверхности Земли, неровного по толщине и неоднородного по латерали верхнего слоя, т.е. когда ни по форме снятых сигналов, ни по результатам Ш-инверсий данных целевой объект выявить невозможно.
ЛИТЕРАТУРА
1. Kaminski V., Prikhodko A., Oldenburg D. Using ERA low frequency E-field profiling and UBC 3D frequency-domain inversion to delineate and discover a mineralized zone in Porcupine district, Ontario, Canada // SEG Technical Program Expanded Abstracts. - 2011. -Vol. 30 (1). - P. 1262-1266. - doi: 10.1190/1.3627433.
2. Parametric 3D inversion of airborne time domain electromagnetics / M.S. McMillan, D.W. Oldenburg, E. Haber, C. Schwarzbach // ASEG Extended Abstracts. - 2015. -Vol. 2015 (1). - P. 1-5. - doi: 10.1071/ASEG2015ab101.
3. Haber E., Schwarzbach C. Parallel inversion of large-scale airborne time-domain electromagnetic data with multiple OcTree meshes // Inverse Problems. - 2014. - Vol. 30, N 5. - P. 055011. - doi: 0266-5611/14/055011.
4. Liu Y., Yin C. 3D inversion for multipulse airborne transient electromagnetic data // Geophysics. - 2015. - Vol. 81 (6). - P. E401-E408. - doi: 10.1190/geo2015-0481.1.
5. Oldenburg D.W., Haber E., Shekhtman R. Three dimensional inversion of multisource time domain electromagnetic data // Geophysics. - 2013. - Vol. 78 (1). - P. E47-E57. -doi: 10.1190/geo2012-0131.1.
6. Yang D., Oldenburg D.W. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit // Geophysics. - 2012. -Vol. 77, N 2. - P. B23-B34. - doi: 10.1190/geo2011-0194.1.
7. Yang D., Oldenburg D.W., Haber E. 3-D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings // Geophysical Journal International. - 2014. - Vol. 196, N 3. - P. 1492-1507. - doi: 10.1093/gji/ggt465.
8. Dehiya R.S.A., Gupta P.K., Israil M. Interpretation of CSEM data using 2D block inversion algorithm // Extended Abstract, 22nd EM Induction Workshop. - Weimar, Germany, 2014. - P. 4.
9. Multiple body parametric inversion of frequency- and time-domain airborne electromagnetics / M.S. McMillan, C. Schwarzbach, E. Haber, D.W. Oldenburg // SEG Technical Program Expanded Abstracts. - 2016. - Vol. 35 - P. 846-851. -doi: 10.1190/segam2016-13868448.1.
10. 3D parametric hybrid inversion of time-domain airborne electromagnetic data / M.S. McMillan, C. Schwarzbach, E. Haber, D. W. Oldenburg // Geophysics. - 2015. -Vol. 80, N 6. - P. K25-K36. - doi: 10.1190/geo2015-0141.1.
11. Application of the marine circular electric dipole method in high latitude Arctic regions using drifting ice floes / V. Mogilatov, M. Goldman, M. Persova, Yu. Soloveichik, Yu. Koshkina, O. Trubacheva, A. Zlobinskiy // Journal of Applied Geophysics. - 2016. - Vol. 135. - P. 1731. - doi: 10.1016/j.jappgeo.2016.08.007.
12. Geometrical nonlinear 3D inversion of airborne time domain EM data / M.G. Persova, Yu.G. Soloveichik, Yu.I. Koshkina, D.V. Vagin, O.S. Trubacheva // Near Surface Geoscience 2016 - First Conference on Geophysics for Mineral Exploration and Mining. Extended abstract. - Barcelona, 2016. - doi: 10.3997/2214-4609.201602114.
13. Методы и алгоритмы восстановления трехмерной структуры проводимости и поляризуемости среды по данным электромагнитных зондирований на основе ко-нечноэлементного 3D-моделирования / М.Г. Персова, Ю.Г. Соловейчик, Г.М. Тригу-бович, М.Г. Токарева // Физика Земли. - 2013. - Т. 3. - С. 30-45.
14. Singh A.D.R., Gupta P.K., Israil M. Development of block Inversion algorithm and its comparison with cell inversion schemes // Extended Abstract, 22nd EM Induction Workshop. -Weimar, Germany, 2014. - P. 4.
15. The topography effect on the airborne EM data / M.G. Persova, Yu.G. Soloveichik, D.V. Vagin, D.S. Kiselev, Yu.I. Koshkina, I.I. Patrushev, E.I. Simon // Saint Petersburg 2018. Innovations in Geosciences - Time for Breakthrough, St. Petersburg, Russia, 9-12 April 2018. - St. Petersburg, 2018. - doi: 10.3997/2214-4609.201800314.
16. 3D time-domain airborne EM forward modeling with topography / C. Yin, Y. Qi, Y. Liu, J. Cai // Journal of Applied Geophysics. - 2016. - Vol. 134. - P. 11-22. -doi: 10.1016/j.jappgeo.2016.08.002.
17. Multidimensional processing of the airborne EM data in the complex media / M.G. Persova, Yu.G. Soloveichik, D.V. Vagin, P.A. Domnikov, D.S. Kiselev, Yu.I. Koshkina, E.I. Simon // Engineering and mining geophysics 2018: 14 conference and exhibition, Kazakhstan, Almaty, 23-27 Apr. 2018. - Almaty, 2018. - Art. У08-01 (46459). - doi: 10.3997/22144609.201800541.
18. Соловейчик Ю.Г., Рояк М.Э., Персова М.Г. Метод конечных элементов для решения скалярных и векторных задач. - Новосибирск: Изд-во НГТУ, 2007. - 896 с.
19. Применение неконформных сеток с шестигранными ячейками для 3D-моделирования технологий аэроэлектроразведки / М.Г. Персова, Ю.Г. Соловейчик, Д.В. Вагин,
Д.С. Киселев, Н.В. Кондратьев, Ю.И. Кошкина, О.С. Трубачева // Доклады Академии наук высшей школы Российской Федерации. - 2018. - № 1 (38). - C. 64-79. -doi: 10.17212/1727-2769-2018-1-64-79. 20. Finite-element solution to multidimensional multisource electromagnetic problems in the frequency domain using non-conforming meshes / Yu.G. Soloveichik, M.G. Persova, P.A. Domnikov, Yu.I. Koshkina, D.V. Vagin // Geophysical Journal International. - 2018. -Vol. 212, N 3. - P. 2159-2193. - doi: 10.1093/gji/ggx530.
PARAMETERIZATION OF A GEOELECTRICAL MODEL IN AIRBORNE ELECTROMAGNETIC PROBLEMS WITH APPLICATIONS TO COMPLEX MEDIA INCLUDING TOPOGRAPHY AND VARIED THICKNESS LAYERS
Kiselev D.S., Kondratyev N.V., Koshkina Yu.I., Vagin D.V., Persova M.G., Soloveichik Yu.G.
Novosibirsk State Technical University, Novosibirsk, Russia
The paper describes the method of the geometric 3-D inversion of airborne electromagnetic time-domain data in significantly inhomogeneous geological media which are characterized by sharp changes in the Earth relief, varying thickness of layers with different conductivity, and laterally inhomogeneous conductivity in the layers overlapping local target objects. In these situations it is impossible to "visually" (using only signal shapes) determine the presence and location of the target object. In addition, neglecting the complex shape of the Earth relief and borders between other layers of the geoelectrical model can either lead to a failure in finding target objects or result in erroneous anomalies. The method is based on the use of a special parameterization which includes physical properties of structural parts of the geoelectrical model and their geometric characteristics. The main feature is the parameterization of the curved borders between the geoelectrical model layers with the use of control points defining the smooth deformation of these borders during the nonlinear 3-D inversion process. In order to represent these surfaces, the bicubic splines are used. The workability of the proposed approach is shown by the example of processing synthetic data obtained for a complex geoelectrical model with the topography and a laterally inhomogeneous overlapping layer with varied thickness. This model has been constructed as generalization on the basis of the results of interpretations of the electromagnetic sounding data obtained on several areas in Eastern Siberia. It is shown that after several iterations, the location (both in the plan and in the depth) of the target object in the form of a kimberlite pipe and its main characteristics have been determined with quite high accuracy.
Keywords: airborne electromagnetic survey, 3-D inversion, topography, electromagnetic field. DOI: 10.17212/1727-2769-2018-4-77-92
REFERENCES
1. Kaminski V., Prikhodko A., Oldenburg D. Using ERA low frequency E-field profiling and UBC 3D frequency-domain inversion to delineate and discover a mineralized zone in Porcupine district, Ontario, Canada. SEG Technical Program Expanded Abstracts, 2011, vol. 30 (1), pp. 1262-1266. doi: 10.1190/1.3627433.
2. McMillan M.S., Oldenburg D.W., Haber E., Schwarzbach C. Parametric 3D inversion of airborne time domain electromagnetics. ASEG Extended Abstracts, 2015, vol. 2015 (1), pp. 15. doi: 10.1071/ASEG2015ab101.
3. Haber E., Schwarzbach C. Parallel inversion of large-scale airborne time-domain electromagnetic data with multiple OcTree meshes. Inverse Problems, 2014, vol. 30, no. 5, p. 055011. doi: 0266-5611/14/055011.
4. Liu Y., Yin C. 3D inversion for multipulse airborne transient electromagnetic data. Geophysics, 2015, vol. 81 (6), pp. E401-E408. doi: 10.1190/geo2015-0481.1.
5. Oldenburg D.W., Haber E., Shekhtman R. Three dimensional inversion of multisource time domain electromagnetic data. Geophysics, 2013, vol. 78 (1), pp. E47-E57. doi: 10.1190/ geo2012-0131.1.
6. Yang D., Oldenburg D.W. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit. Geophysics, 2012, vol. 77, no. 2, pp. B23-B34. doi: 10.1190/geo2011-0194.1.
7. Yang D., Oldenburg D.W., Haber E. 3-D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings. Geophysical Journal International, 2014, vol. 196, no. 3, pp. 1492-1507. doi: 10.1093/gji/ggt465.
8. Dehiya R.S.A., Gupta P.K., Israil M. Interpretation of CSEM data using 2D block inversion algorithm. Extended Abstract, 22nd EM Induction Workshop, Weimar, Germany, 2014, p. 4.
9. McMillan M.S., Schwarzbach C., Haber E., Oldenburg D.W. Multiple body parametric inversion of frequency- and time-domain airborne electromagnetics. SEG Technical Program Expanded Abstracts, 2016, vol. 35, pp. 846-851. doi: 10.1190/segam2016-13868448.1.
10. McMillan M.S., Schwarzbach C., Haber E., Oldenburg D.W. 3D parametric hybrid inversion of time-domain airborne electromagnetic data. Geophysics, 2015, vol. 80, no. 6, pp. K25-K36. doi: 10.1190/geo2015-0141.1.
11. Mogilatov V., Goldman M., Persova M., Soloveichik Yu., Koshkina Yu., Trubacheva O., Zlobinskiy A. Application of the marine circular electric dipole method in high latitude Arctic regions using drifting ice floes. Journal of Applied Geophysics, 2016, vol. 135, pp. 17-31. doi: 10.1016/j.jappgeo.2016.08.007.
12. Persova M.G., Soloveichik Yy.G., Koshkina Yu.I., Vagin D.V., Trubacheva O.S. Geometrical nonlinear 3D inversion of airborne time domain em data. Near Surface Geoscience 2016 - First Conference on Geophysics for Mineral Exploration and Mining. Extended abstract, Barcelona, 2016. doi: 10.3997/2214-4609.201602114.
13. Persova M.G., Soloveichik Yu.G., Trigubovich G.M., Tokareva M.G. Methods and algorithms for reconstructing three-dimensional distributions of electric conductivity and polarization in the medium by finite-element 3D modeling using the data of electromagnetic sounding. Izvestiya, Physics of the Solid Earth, 2013, vol. 49 (3), pp. 329-343. doi: 10.1134/S1069351313030117. Translated from Fizika Zemli, 2013, vol. 3, pp. 30-45.
14. Singh A.D.R., Gupta P.K., Israil M. Development of block Inversion algorithm and its comparison with cell inversion schemes. Extended Abstract, 22nd EM Induction Workshop, Weimar, Germany, 2014, p. 4.
15. Persova M.G., Soloveichik Yu.G., Vagin D.V., Kiselev D.S., Koshkina Yu.I., Patrushev I.I., Simon E.I. The topography effect on the airborne EM data. Saint Petersburg 2018. Innovations in Geosciences — Time for Breakthrough, St. Petersburg, Russia, 9-12 April 2018, doi: 10.3997/2214-4609.201800314.
16. Yin C., Qi Y., Liu Y., Cai J. 3D time-domain airborne EM forward modeling with topography. Journal of Applied Geophysics, 2016, vol. 134, pp. 11-22. doi: 10.1016/j.jappgeo.2016.08.002.
17. Persova M.G., Soloveichik Yu.G., Vagin D.V., Domnikov P.A., Kiselev D.S., Koshkina Yu.I., Simon E.I. Multidimensional processing of the airborne EM data in the complex media. Engineering and mining geophysics 2018: 14 conference and exhibition, Kazakhstan, Almaty, 23-27 Apr. 2018, art. У08-01 (46459). doi: 10.3997/22144609.201800541.
18. Soloveichik Yu.G., Royak M.E., Persova M.G. Metod konechnykh elementov dlya resheniya skalyarnykh i vektornykh zadach [The finite element method for the solution of scalar and vector problems]. Novosibirsk, NSTU Publ., 2007. 896 p.
19. Persova M.G., Soloveichik Y.G., Vagin D.V., Kiselev D.S., Kondratyev N.V., Kosh-kina Y.I., Trubacheva O.S. Primenenie nekonformnykh setok s shestigrannymi yacheikami dlya 3D-modelirovaniya tekhnologii aeroelektrorazvedki [Application of non-conforming meshes with hexahedral cells for3D modelling of airborne electromagnetic technologies]. Doklady Akademii nauk vysshei shkoly Rossiiskoi Federatsii - Proceedings of the Russian higher school Academy of sciences, 2018, no. 1 (38), pp. 64-79. doi: 10.17212/1727-27692018-1-64-79.
20. Soloveichik Yu.G., Persova M.G., Domnikov P.A., Koshkina Yu.I., Vagin D.V. Finite-element solution to multidimensional multisource electromagnetic problems in the frequency domain using non-conforming meshes. Geophysical Journal International, 2018, vol. 212, no. 3, pp. 2159-2193. doi: 10.1093/gji/ggx530.
СВЕДЕНИЯ ОБ АВТОРАХ
Киселев Дмитрий Сергеевич - родился в 1990 году, аспирант, ассистент кафедры прикладной математики Новосибирского государственного технического университета. Область научных интересов: конечноэлементное моделирование. Опубликовано 20 научных работ. (Адрес: 630073, Россия, Новосибирск, пр. Карла Маркса, 20. E-mail: harlequin_00@mail.ru).
Kiselev Dmitry Sergeevich (b. 1990) - a postgraduate student, an assistant lecturer at the department of applied mathematics, Novosibirsk State Technical University. His research interests are currently focused on finite element modeling. He is the author of 20 scientific papers. (Address: 20, Karl Marx Av., Novosibirsk, 630073, Russia. E-mail: harlequin_00@mail.ru).
Кондратьев Николай Владимирович - родился в 1991 году, ассистент кафедры прикладной математики Новосибирского государственного технического университета. Область научных интересов: разработка и оптимизация параллельных программ решения СЛАУ, полученных в результате конечноэлементной аппроксимации в задачах электромагнетизма в системах с общей и разделяемой памятью. Опубликовано 7 научных работ. (Адрес: 630073, Россия, Новосибирск, пр. Карла Маркса, 20. E-mail: kondratyev@ami.nstu.ru).
Kondratyev Nikolay Vladimirovich (b. 1991) - an assistant lecturer at the department of applied mathematics, Novosibirsk State Technical University. His research interests are currently focused on the development and optimizations of parallel programs, acceleration of numerical modeling on CPU and GPU devices, development and implementation of distributed calculation systems for numerical modeling. He is the author of 7 scientific papers. (Address: 20, Karl Marx Av., Novosibirsk, 630073, Russia. E-mail: kondratyev@ami.nstu.ru).
Кошкина Юлия Игоревна - родилась в 1990 году, канд. техн. наук, доцент кафедры прикладной математики Новосибирского государственного технического университета. Область научных интересов: разработка и реализация алгоритмов интерпретации данных электромагнитных зондирований. Опубликовано 30 научных работ. (Адрес: 630073, Россия, Новосибирск, пр. Карла Маркса, 20. E-mail: koshkina_yui@mail.ru).
Koshkina Yulia Igorevna (b. 1990) - Candidate of Sciences (Eng.), an associate professor at the department of applied mathematics, Novosibirsk State Technical University. Her research interests are currently focused on the development and implementation of algorithms for electromagnetic sounding data interpretation. She is the author of 30 scientific papers. (Address: 20, Karl Marx Av., Novosibirsk, 630073, Russia. E-mail: koshkina_yui@mail.ru).
Вагин Денис Владимирович - родился в 1985 году, канд. техн. наук, доцент кафедры прикладной математики Новосибирского государственного технического университета. Область научных интересов: конечно-элементное моделирование электромагнитных полей. Опубликовано более 70 научных работ. (Адрес: 630073, Россия, Новосибирск, пр. Карла Маркса, 20. E-mail: vdv_wk@mail.ru).
Vagin Denis Vladimirovich (b. 1985) - Candidate of Sciences (Eng.), an associate professor at the department of applied mathematics, Novosibirsk State Technical University. His research interests are currently focused on finite element modeling of electromagnetic fields. He is the author of more than 70 scientific papers. (Address: 20, Karl Marx Av., Novosibirsk, 630073, Russia. Email: vdv_wk@mail.ru).
Персова Марина Геннадьевна - родилась в 1978 году, д-р техн. наук, профессор, профессор кафедры прикладной математики Новосибирского государственного технического университета. Область научных интересов: конечноэлементное моделирование и наукоемкое программное обеспечение для моделирования и сопровождения технологий. Опубликовано более 150 научных работ. (Адрес: 630073, Россия, Новосибирск, пр. Карла Маркса, 20. E-mail: persova@ ami.nstu.ru).
Persova Marina Gennad'evna (b. 1978) - Doctor of Sciences (Eng.), professor, professor at the applied mathematics department, Novosibirsk State Technical University. Her research interests are currently focused on finite element modeling and knowledge-intensive software for modeling and supporting technologies. She is the author of more than 150 scientific papers. (Address: 20, Karl Marx Av., Novosibirsk, 630073, Russia. E-mail: persova@ ami.nstu.ru).
Соловейчик Юрий Григорьевич - родился в 1957 году, д-р техн. наук, профессор, член-корреспондент САН ВШ, заведующий кафедрой прикладной математики Новосибирского государственного технического университета. Область научных интересов: конечноэлементное моделирование и наукоемкое программное обеспечение для моделирования и сопровождения технологий. Опубликовано более 150 научных работ. (Адрес: 630073, Россия, Новосибирск, пр. Карла Маркса, 20. E-mail: soloveychik@ami.nstu.ru).
Soloveichik Yuri Grigorievich (b. 1957) - Doctor of Sciences (Eng.), professor, head of the department of applied mathematics, Novosibirsk State Technical University. His research interests are currently focused on finite element modeling and knowledge-intensive software for modeling and supporting technologies. He is the author of more than 150 scientific papers. (Address: 20, Karl Marx Av., Novosibirsk, 630073, Russia. E-mail: soloveychik@ami.nstu.ru).
Статья поступила 25 октября 2018 г.
Received October 25, 2018
ft
To Reference:
Kiselev D.S., Kondratyev N.V., Koshkina Yu.I., Vagin D.V., Persova M.G., Soloveichik Yu.G.
0 parametrizatsii geoelektricheskoi modeli v zadachakh aeroelektrorazvedki v sredakh s rel'efom
1 sloyami peremennoi tolshchiny [Parameterization of a geoelectrical model in airborne electromagnetic problems with applications to complex media including topography and varied thickness layers]. Doklady Akademii nauk vysshei shkoly Rossiiskoi Federatsii - Proceedings of the Russian higher school Academy of sciences, 2018, no.4 (41), pp. 77-92. doi: 10.17212/17272769-2018-4-77-92.
