ФИЗИКА ЗЕМЛИ, 2013, № 3, с. 30-45
УДК 550.837:517.958
МЕТОДЫ И АЛГОРИТМЫ ВОССТАНОВЛЕНИЯ ТРЕХМЕРНОЙ СТРУКТУРЫ ПРОВОДИМОСТИ И ПОЛЯРИЗУЕМОСТИ СРЕДЫ ПО ДАННЫМ ЭЛЕКТРОМАГНИТНЫХ ЗОНДИРОВАНИЙ НА ОСНОВЕ КОНЕЧНОЭЛЕМЕНТНОГО 3D-МОДЕЛИРОВАНИЯ © 2013 г. М. Г. Персова1, Ю. Г. Соловейчик1, Г. М. Тригубович2, М. Г. Токарева1
Новосибирский государственный технический университет, г. Новосибирск 2Сибирский научно-исследовательский институт геологии, геофизики и минерального сырья, г. Новосибирск
Поступила в редакцию 21.05.2012 г.
В работе рассматриваются методы и алгоритмы автоматической ЗБ-инверсии данных становления поля и вызванной поляризации, основанные на использовании конечноэлементного ЗБ-модели-рования. На синтетических данных, полученных с использованием ЗБ-моделирования, приводятся примеры работы алгоритмов автоматической ЗБ-инверсии.
БО1: 10.7868/80002333713030113
ВВЕДЕНИЕ
Использование трехмерных подходов при интерпретации электроразведочных данных позволяет существенно повысить качество восстановления характеристик поисковых объектов. Одна из главных проблем интерпретации практических данных — появление ложных аномалий и пропуск реальных, что довольно часто случается при использовании упрощенных методов интерпретации из-за неучета или неправильного учета неоднородностей верхней части разреза, перекрывающих поисковые объекты. Соответствующие примеры приведены, например, в работе [Тригубович и др., 2009]. В этой работе также приведены результаты ЗБ-инверсий практических данных зондирований становлением поля, полученных на некоторых площадях Восточной Сибири и Дальнего Востока, которые были выполнены с использованием программного комплекса ОеоБМ [Персова и др., 2011а]. Этот программный комплекс, позволяющий выполнять ЗБ-моделирование геоэлектромагнитных полей для различных источников, основан на конечно-элементном моделировании по технологии с выделением поля [Соловейчик и др., 1998; Персова и др., 2011а] и использовании оптимизированных автоматически строящихся быстро разрежающихся нерегулярных параллелепипеидальных сеток [Соловейчик и др., 2007].
Однако эти инверсии выполнялись "вручную", при построении трехмерных моделей среды по некоторым площадям операторам требовалось решать до нескольких тысяч трехмерных
задач. Безусловно, эта работа очень трудоемка и требует высокой квалификации операторов. Поэтому разработка алгоритмов автоматических 3Э-инверсий является актуальной проблемой, которая в настоящее время довольно широко обсуждается в научных публикациях. Описанные в работах [Жданов, 2007; Farquharson er al., 2002; Gribenko, Zhdanov, 2007; Haber et al., 2007; Mackie et al., 2001; Newman, Boggs, 2004; Newman, Com-mer, 2005; Siripunvaraporn et al., 2005; Sasaki, 2001; Sasaki et al., 2010; Varentsov, 2005; Yoshioka, Zhdanov, 2005; и др.] подходы к решению многомерных обратных задач основаны на том, что исследуемая область представляется в виде набора пробных объектов, т.е. в виде некоторой сетки, в каждой ячейке которой (пробном объекте) ищется значение удельной проводимости. Этот подход при выполнении полной нелинейной инверсии (т.е. без линеаризации зависимости поля от проводимости трехмерных объектов как при вычислении теоретических значений сигналов для очередного приближения, так и при расчете производных отклонений по параметру) является очень вычислительно затратным. При этом одним из критериев качества такого рода подходов является возможность использования сеток с достаточно мелкими пробными объектами (с целью более точного описания границ локальных неоднород-ностей). В данной работе будет рассмотрен несколько модифицированный подход к решению обратной задачи, основанный на двухэтапной процедуре, на первом шаге которой будет искаться стартовая геоэлектрическая модель распределения проводимости с относительно крупными
пробными объектами, а на втором шаге будут уточняться параметры локальных неоднородно-стей, сформированных из пробных объектов с близкими значениями удельной проводимости.
Необходимо также отметить, что в подавляющем большинстве работ, посвященных разработке методов выполнения многомерных инверсий, рассматриваются задачи, в которых электромагнитное поле возбуждается гармоническим током, в том числе это задачи аэроэлектроразведки, задачи магнитотеллурики, индукционного каротажа. И довольно мало появляется работ, посвященных методам выполнения трехмерных инверсий во временной области. Кроме того, несмотря на то, что в последних работах для расчета теоретических сигналов на очередной итерации уже довольно часто используется точное моделирование, вычисление производных по параметрам чаще всего выполняется с использованием определенных упрощений, что может негативно влиять как на сходимость, так и на результаты инверсии.
В данной работе мы рассмотрим подходы к 3Э-инверсии данных электромагнитных зондирований во временной области, полученных методом становления поля и методом вызванной поляризации, в которых решение прямой задачи как при расчете теоретических сигналов для очередного приближения, так и при расчете производных по параметрам будет основано на точном моделировании с использованием метода конечных элементов (МКЭ).
МАТЕМАТИЧЕСКИЕ МОДЕЛИ ДЛЯ РЕШЕНИЯ ПРЯМЫХ ЗАДАЧ МОДЕЛИРОВАНИЯ ПРОЦЕССОВ СТАНОВЛЕНИЯ ПОЛЯ И ВЫЗВАННОЙ ПОЛЯРИЗАЦИИ И АЛГОРИТМЫ 3Э-ИНВЕРСИИ
Вначале рассмотрим кратко математический аппарат, на котором основаны предлагаемые подходы к 3Э-инверсии данных становлением поля. Параметры трехмерной геоэлектрической модели будут определяться на основании минимизации функционала
ь к
Фа (Ь) = ХЕК86« (Ь))2 + С(Ь - Ь)2,
(1)
1=1 к=1
гуляризирующую добавку; С — матрица, на главной диагонали которой расположены параметры регуляризации (и все внедиагональные элементы С равны 0); ю!к — некоторые веса, отражающие уровень погрешности при приеме сигнала в 1-м приемнике и масштаб изменения принимаемого сигнала по времени. В качестве весов ю!к могут быть взяты, например, значения, обратные к значениям ЭДС в центральной точке генераторной петли. Вектор искомых параметров Ь (размерностью М) включает в себя значения удельной проводимости искомых объектов (заданных в исследуемой области в виде прямоугольных параллелепипедов), а также может включать границы этих объектов по х, у и г.
Линеаризуем отклонения 5б;,
м
т=1
8е№ (Ь) ~ 8е,к (Ь0) + Х^ АЬт,
т
(2)
где Ь0 — вектор параметров, полученный на преды,, . ,о д (5е 1к) дущей итерации, АЬт = Ьт — Ьт, а —-— произ-
дЬт
водные, отражающие влияние изменения т-го параметра в 1-й приемной петле в момент времени с номером к, и подставим представление (2) в функционал (1):
Фа (Ь) =
ь к
1=1 к=1
м
ХХ Ю »5е,к (Ь°) + ЮкХ
д(8Ек)
т=1
дЬт,
ЛЬ.
(3)
в(ь0 - Ь + ЛЬ)2
В результате минимизации функционала (3) по ЫЬт получаем СЛАУ вида
(А + в)АЬ = I - в (Ь0 - Ь),
(4)
где I — элементы матрицы А и вектора правой части Г определяются соотношениями
ь к
А = ХХ((°1к)
1=1 к=1
\ 2 5 (5 к )д (5е 1к) дЬ1 дЬ.. '
ь к
где Ь — вектор искомых параметров; Ъг1к = § 1к — е?к, &1к — сигналы ЭДС, зарегистрированные в 1-м приемнике в моменты времени 1к, § 1к — теоретические сигналы, полученные в результате решения
прямой трехмерной задачи; Ь — вектор некоторых фиксированных параметров, определяющих ре-
^ =-ХХК )2 8% (Ь0 , 1 = 1... м.
1=1 к=1 1
Решение прямой задачи при вычислении б 1к (Ь) выполняется с использованием векторного МКЭ для математической модели, основанной на так называемой технологии многоэтапного выделения поля [Соловейчик и др., 2007; Персова и др., 2011в]. При этом прямая сумма получаемых полей равна искомому полю полной трехмерной задачи. Алгоритм этого метода заключается в следу-
ющем. На первом шаге вычисляется поле горизонтально-слоистой среды (как решение задачи меньшей размерности). Затем на втором шаге вычисляется поле влияния первого объекта, т.е. решается трехмерная задача в области, представляющей собой горизонтально-слоистую среду с одной трехмерной неоднородностью. На третьем шаге вычисляется поле влияния второго объекта относительно среды, содержащей первый объект, т.е. решается задача в среде, содержащей два трехмерных объекта, но при этом только второй объект является источником аномального поля. Такая процедура повторяется для всех ЗЭ-объектов модели, и на последнем шаге выполняется расчет поля влияния последнего объекта относительно среды, содержащей все остальные трехмерные объекты. Математическая модель для расчета поля влияния трехмерного объекта имеет вид [Персова и др., 2011б]
rot
— rotÄ'
vMq J
+ a
sä! dt
= |a
I 3D Q\p
la-a - )E
3D Q
(5)
Значения
также ищутся с использова-
3Д 0 _0
где а - и Е — распределения проводимости и напряженности электрического поля в трехмерной среде, поле для которой было рассчитано на предыдущем этапе алгоритма и относительно которой на данном этапе вычисляется поле влияния очередного объекта. Заметим, что в уравнении (5) а — распределение проводимости в трехмерной среде, содержащей трехмерные объекты, поля влияния которых вычислялись на предыдущих этапах алгоритма, и объекта, поле влияния которого вычисляется на текущем этапе. Поэтому а Ф азв_0 только в месте расположения текущего объекта. На последнем этапе работы алгоритма а — распределение проводимости в трехмерной среде, соответствующей полной трехмерной модели.
Аномальная составляющая напряженности электрического поля на каждом шаге работы алго-
—*а
а дА
ритма определяется в виде Е =--. На втором
дt
же шаге работы алгоритма, когда рассчитывается
3Ю 0
поле влияния первого объекта, в качестве а - и
3Б_ 0 Ы Ы
Е берутся а и Е — распределения проводимости и напряженности электрического поля во вмещающей горизонтально-слоистой среде. Математические модели для расчета Е для различных источников электромагнитного поля приведены в работе [Персова и др., 2011а]. Для аппроксимации по времени уравнения (5) используется трехслойная неявная схема с увеличивающимся шагом по времени, а для аппроксимации по пространству — векторные базисные функции перво-
го порядка (edge-элементы) [Bossavit, 1998; Соловейчик и др., 2007
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.