CC BY
9Научный вестник НГТУ. - 2013. -№ 2(51)
УДК 519.63
Алгоритмы 31)-инверсии данных зондирований становлением поля с использованием борновских приближений
М.Г. ПЕРСОВА, Е.И. СИМОН, Ю.Г. СОЛОВЕЙЧИК, Ю.И. КОШКИНА
Рассматриваются алгоритмы ЗО-инверсии данных зондирований становлением поля, базирующиеся на решении линеаризованной обратной задачи на основе борновских приближений. Приводится математический аппарат и результаты работы этих алгоритмов на синтетических аналогах полевых данных, полученных с использованием математического ЗО-моделирования.
Ключевые слова: ЗО-инверсии, борновские приближения, ЗО-моделирование. ВВЕДЕНИЕ
Создание программного обеспечения, реализующего ЗБ-инверсии данных электромагнитных зондирований, является на сегодняшний день самой актуальной проблемой геоэлектрики 11 —41. Стандартным подходом при реализации алгоритмов ЗБ-инверсии данных электромагнитных зондирований является разбиение изучаемого объема среды на ячейки и поиск коэффициента удельной электрической проводимости в каждой из ячеек. Очевидно, что главным препятствием развития и внедрения этих алгоритмов в практику электроразведочных исследований является их вычислительная трудоемкость, и поэтому чаще всего при реализации этих алгоритмов используют упрощенные математические модели, что может негативно влиять как на сходимость, так и на результаты инверсии.
Кроме того, повышения качества результатов ЗБ-инверсии, как правило, пытаются добиться путем использования более мелких разбиений, что довольно часто помимо резкого роста вычислительных затрат приводит к нахождению одного из «эквивалентных» решений с довольно «пестрой» картиной распределения удельного сопротивления, что не только не облегчает, а наоборот, затрудняет определение морфологии и удельного сопротивления целевого объекта.
Решением данной проблемы может являться, например, использование многоэтапных процедур, в которых после первого этапа, где с помощью ячеистой структуры ищется стартовое распределение параметров, выполняется второй этап, на котором уточняются параметры локальных неоднородностей, сформированных по результатам инверсии в ячеистой структуре с достаточно крупными ячейками. Кроме того, на первых этапах могут использоваться различные варианты быстрых инверсий для получения стартовых геоэлектрических моделей, уточнение которых на последующих этапах будет осуществляться на основе нелинейных ЗБ-инверсий с использованием точных (без упрощений) моделей геоэлектромагнитного поля.
Данная работа посвящена разработке такого рода процедур выполнения ЗБ-инверсий, которые могут быть использованы на этапе построения стартовых геоэлектрических моделей среды, представляющей собой набор 3 Г) - м а крообъе кто в, уточнение параметров которых вы-
Статпья получена 4 декабря 2012 г.
Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации.
полняется в ходе последующего применения процедур нелинейных ЗБ-инверсий, а также для уточнения количества и местоположения макрообъектов на промежуточных итерациях нелинейной 3 D-ипвсрсии.
1. МАТЕМАТИЧЕСКИЙ АППАРАТ ЗО-ИНВЕРСИЙ ДАННЫХ ЗОНДИРОВАНИЙ СТАНОВЛЕНИЕМ ПОЛЯ С ИСПОЛЬЗОВАНИЕМ БОРНОВКИХ ПРИБЛИЖЕНИЙ
Электромагнитное нестационарное поле в однородной по магнитной проницаемости среде, вызванное круговой генераторной петлей с током, находящейся в плоскости z — const, может быть полностью описано решением следующей начально-краевой задачи:
—rot (rot 6 ) + а— = Jст ,в X, 1\ 2 О,
M0 V 'dt Ir '
(1)
"гСТ
где а - электрическая проводимость, |Д0- магнитная проницаемость вакуума, .7 - плотность стороннего тока (определяемая током в генераторной петле), Г - граница расчетной области X .В качестве начального условия 6 берется решение стационарной задачи
It 2t0
—rot (rot 6 ) = Jст ,в X, 6 2 0, M0 lr
(2)
если в качестве источника задана функция Хевисайда (на практике это означает достаточно длинный импульс, то есть такой, что к моменту выключения тока в генераторной петле электромагнитное поле уже полностью установилось).
При использовании технологии разделения поля на нормальную и аномальную состав— — —0 —а
ляющие вектор А представляется в виде суммы: А — А + А . При численном моделировании геоэлектромагнитных полей в качестве нормального поля используется поле источника в горизонтально-слоистой среде [5-8], которое можно получить через решение двумерной задачи меньшей размерности. Для нахождения поля влияния трехмерных неоднородностей необходимо решать трехмерное векторное уравнение
—rot (rot 6а ) + а— = (ст-ст° )@ °,в X, 6
M o
2 0,
Е o 36 где Е 2 -
dt
6а
2 0.
12t0
Аномальная составляющая Аа вектор-потенциала А так>ке удовлетворяет уравнению
J_ M0
rot
(rot 6а )-
а
36а dt
: (а — а0 )Е ,
£ 8А
где Ь —---напряженность суммарного электрического поля.
дt
В соответствии с методом интегральных уравнений, решение уравнения (4) в любой точке пространства может быть получено как значение реакции горизонтально-слоистой среды в
этой точке на поле, возникающее в аномалии. Интегральное представление для вектора Аа
6а (',у,t)2 J JGу, z,' ',у',z' , t — т))
—даП
\ dE
а — а))—('у ,zT)dx'dy'dz'dт , (5)
где G(x,y,z,x',y',z',t — т) - тензор Грина, компоненты которого определяются как реакция
электромагнитного поля (в виде значений вектор-потенциала Я") горизонтально-слоистой среды с удельной проводимостью ст° в момент времени т в точке с координатами {х,у, z j на соответствующим образом сориентированный точечный источник (электрический диполь), расположенный в точке (x',y',z'), действовавший бесконечно долго до момента времени
t = 0 и выключенный в этот момент (7 = 0).
Заменим в соотношении (5) стоящее под интегралом поле Е на поле Еп среды с проводимостью а", в результате получим линеаризованное поле
q t . ,
Аа (х,у,z,t)= J J G(x,у,z,x',у',z',t -t)((j" -ak J-(xу',z',x)dx'dy'dz' dx, (6)
-»Xt dx
где Q/( - подобласти, в которых проводимость а, равная а/(, отличается от проводимости ст" среды, полученной на п -ой итерации процедуры восстановления проводимости изучаемой трехмерной среды. Поле Еп может быть найдено из решения соответствующей начально-краевой задачи, аналогичной задаче для уравнения (3).
Нестационарное электромагнитное поле, вызванное горизонтальным электрическим диполем, согласно подходу, предложенному в работах B.C. Могилатова [9], описывается уравнениями
92 х dx п --— + У2 х + ца-2 0,
9z2 dt
1 9V ) .2V 9V n
--1 + У — + ц-2 0,
dz 1а dz ) а dt
(&) (8)
с краевыми условиями
[цх d - 2 „ 2о, [к ]
|-2Xq (t), i 21,
0,
i ФI,
X ^ 0,
[V ]
f-2q,(')'i 2 [V/а]
0, i Ф У,
z2Zi 2 0, V ^ 0, z
где - координаты, соответствующие границам слоев, - координата, соответствующая положению источника.
Компоненты электромагнитного поля II - рассчитываются через функции X и V
Н2 2 Ы'0 —
д ду
-1-JJ0 (Уг)X(y)dУ 4л 0
(9)
где У0(А,г) - функции Бесселя.
Программа, реализующая ЗБ-инверсию данных на основе борновских приближений, полученных по технологии зондирования становлением поля, основана на следующей вы-
числительной схеме. Будем считать, что для нескольких (Р) положений генераторной установки (источника поля) получены сигналы с(7] в + \>2 +... + ¡> = £ приемниках (т.е. для р-то положения источника сигналы с (7] регистрировались в V ¡> приемниках, р — \,...,Р). Обозначим через с//( аномальные сигналы ЭДС, зарегистрированные в 1-м приемнике в момент времени tk [I — \...Ь, к — \... К^. Если для двух положений источника
используется одно и то же положение приемника, то будем считать, что это разные приемники (т.е. они имеют разные номера /). Аномальные теоретические сигналы, полученные в результате решения линеаризованной прямой задачи с использованием борновского приближения, обозначим через .
Исследуемая часть среды разбивается на элементарные подобласти 0.т, т = 1... М, в каждой из которых будут искаться свои значения Дстт =от _ ^ (где ст^ проводимость референтной среды, относительно которой вычисляются аномальные сигналы с//().
Значения с//(, как уже говорилось, ищутся в результате решения линеаризованной прямой задачи в виде:
м
2 X ДатГ 1к >
м=1
(10)
где - значения ЭДС поля, рассчитываемого в референтной среде и порождаемого источниками в виде единичных электрических диполей, заданных в подобластях 0.т с аномальной проводимостью ат .
Значения Аат ищутся в результате минимизации функционала
Фа(Аат ) = !!
121 к 21
' (м &1к
V
Х(т гЛ )
^ т21
ЛЛ
%
//
м
!
"Хат(Д°т) +
т21
м
ХУт X (
(11)
2
т21
1)с1 т
где юц- - некоторые веса, ат и ут - параметры регуляризации, а 1т - множество номеров ячеек, окружающих т -ю ячейку и входящих в ту же подобласть сглаживания. В качестве используются величины, обратные к значениям ЭДС горизонтально-слоистой вмещающей среды. Значения параметров регуляризации «,,„ выбираются максимальными, при которых значение функционала увеличивается не более, чем на 1 %, а также такими, чтобы найденные значения Даш соответствовали положительным значениям <зт . Параметры уш определяются необходимым уровнем гладкости получаемого распределения проводимости а.
Фактически задача минимизации (11) сводится к решению СЛАУ
АД б 2 Г
для вектора Бэ = (Аа1,...,Аат)т
аномальных проводимостеи.
А и вектора Г в (12)
(13)
121 X 21
I К
!
6 2 XX Ц*) гХг/* + «V + X (уV + У^)
(14)
У1 X 21
I К
2ХХ(<%)!%4,}213М.
(15)
У 21 х 21
г1к
£1т ( <5т ) задается равной
^(*)2 )! +()!лх.
(16)
По решению линеаризованной задачи определяются направления, вычисляется новое приближение значений удельной проводимости и рассчитывается функционал невязки нелинейной обратной задачи путем конечноэлементного решения прямой задачи для уравнения,
аналогичного (3). На следующей итерации при поиске очередного направления в качестве Е" и Е" берутся значения напряженности электрического поля, полученные при решении трехмерной задачи для распределения удельной проводимости а" с предыдущей итерации.
По окончании итерационного процесса на основе полученного распределения параметров удельного сопротивления выделяются подобласти ячеек, предположительно соответствующие локальным трехмерным неоднородностям. После этого процесс может быть продолжен с использованием специальной процедуры сглаживания, основанной на близости значений удельного сопротивления внутри выделенных подобластей, путем выбора подходящих значений ут в функционале (11).
Верификация разработанных методов автоматических ЗБ-инверсий выполнялась с использованием различных геоэлектрических моделей, для которых с помощью ЗБ-моделирования были синтезированы аналоги полевых данных. В качестве одного из примеров рассмотрим геоэлектрическую модель, представляющую собой пятислойную вмещающую среду с параметрами \ =230 м, Р! =70 Ом-м, =270 м, р2 =15 Ом-м, /г3 =440 м, Рз =7 Ом-м, /г4 =460 м, р4 =3 Ом-м, р5 =300 Ом-м (/гг- - толщина / -го слоя, рг- - его удельное сопротивление), в которую на разные глубины помещено 3 объекта. Эта модель представлена на рис. 1, где показана квадратная генераторная петля размером 1x1 км2, а также соответствующие приведенному положению генераторной петли приемники (обозначены точками) и траектория перемещения приемно-генераторной установки по трем профилям (на каждом из профилей по четыре положения приемно-генераторной установки).
В результате работы процедуры нелинейной ЗБ-инверсии в двух слоях ячеек было восстановлено удельное электрическое сопротивление, показанное на рис. 2. Всего было сделано 3 итерации, что обеспечило понижение функционала невязки от 24.2 до 2.54. При этом значения ут в функционале (11) принимались равными нулю. Черной пунктирной линией здесь и далее
объектов, расположенных в соответствующих слоях.
р=1 ОМ'М
р=30 Ом-м р=100 Ом-м
Рис. 1. Истинная геоэлектрическая модель
По полученному распределению удельного электрического сопротивления были выделены подобласти, предположительно соответствующие местоположению трехмерных объектов, и процесс ЗБ-инверсии был продолжен, но уже с использованием специальной процедуры сглаживания, обеспечивающей близость значений удельного сопротивления в ячейках, находящихся внутри каждой из выделенных подобластей (т. е. при ут Ф 0 , и при этом соседними считались только ячейки, принадлежащие одной подобласти).
■
Л4
Рис. 2. Распределение удельного электрического сопротивления в слое ячеек от 500 до 940 м (а) и в слое ячеек от 940 до 1400 м (б)
На рис. 3 приведено распределение удельного электрического сопротивления, полученное с использованием специальной процедуры сглаживания.
Рис. 3. Распределение удельного электрического сопротивления в слое ячеек от 500 до 940 м (а) и в слое ячеек от 940 до 1400 м (б), полученные после применения специальной
процедуры сглаживания
Из приведенных на рис. 3 результатов видно, что полученное распределение удельного электрического сопротивления достаточно хорошо соответствует истинной модели.
В качестве следующего примера рассмотрим геоэлектрическую модель, показанную на рис. 4. Эта геоэлектрическая модель имеет более сложную структуру и представляет собой шестислойную вмещающую среду с параметрами \ = 300 м, р| = 300 Ом-м, И2 = 200 м, р2 =5 Ом-м, /?3 =1700 м, р3 =170 Ом-м, /?4 =100 м, р4 =10 Ом-м, /?5 =500 м, р5 = 500 Ом-м, р6 = 5000 Ом-м, в верхнюю часть которой помещено 5 объектов (имитирующих геологические объекты-помехи, с удельными сопротивлениями от 3 до 10 Ом-м), а в нижнюю часть - целевой объект с удельным сопротивлением 2 Ом-м.
а б
Рис. 4. Геоэлектрическая модель, содержащая приповерхностные объекты-помехи и глубинный целевой объект: план (а), объемное изображение (о)
На первом этапе интерпретации данных по комплексу всех данных была подобрана некоторая средняя вмещающая среда с параметрами \ =300 м, Р| =300 Ом-м, /ъ = 200 м, р0 =5.13 Ом-м, /?3 =200 м, р3 =75 Ом-м, /?4 =2850 м, р4 =200 Ом-м, р5 =840 Ом-м. На втором этапе в ячеистой структуре, построенной на уровне верхнего слоя, с использованием разработанной процедуры ЗП-ипвсрсии было подобрано удельное электрическое сопротивление. При выполнении ЗБ-инверсии было сделано 3 итерации, что позволило понизить функционал невязки со значения 55 до 22.7. Соответствующее распределение удельного сопротивления показано на рис. 5, а. На рис. 5, б показано распределение удельного сопротивления после использования процедуры сглаживания проводимости по подобластям.
На следующем этапе подобранные значения удельного электрического сопротивления в ячейках в верхних слоях были зафиксированы, и был осуществлен подбор удельного сопротивления в ячейках, расположенных на глубине.
На рис. 6 приведено распределение удельного электрического сопротивления в нижнем слое с использованием процедуры сглаживания во всей области и в выделенных подобластях.
Для дальнейшего улучшения подбора параметров глубинного объекта требуется выполнение полной нелинейной ЗП-ипвсрсни - процедуры, где направления также ищутся с использованием точного ЗО-моделирования, а не на основе борновских приближений [10].
Тем не менее, разработанная процедура, основанная на борновской инверсии, может дать правильные «подсказки» пользователю для задания стартовой модели в полной нелинейной инверсии - показать места «дефицита» и «избытка» проводимости.
i
-
/
/_.
-40
- 36
- 32
- 28 24 20
- 16 - 12
-IÚÚÚÚ -8ÚÚÚ -6ÚÚÚ -4000 -2000
2000 4000 6000 8000 IÚÚÚÚ -10000 -8000 -6000 -4000 -2000
2000 4000 6000 8000 10000^™ 0.1
-2
I
1000^ 0
Рис. 5. Распределение удельного электрического сопротивления в верхнем слое после трех итераций до использования процедуры сглаживания (а) и после использования процедуры сглаживания в выделенных подобластях (б)
-10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 10000
10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000
Рис. 6. Распределение удельного электрического сопротивления в глубинном слое с использованием процедуры сглаживания во всей области (а) и в выделенных подобластях ( б)
Например, для подобранного в результате описанной выше процедуры ЗО-инверсии распределения удельного сопротивления, можно оценить избыточную (или недостаточную) проводимость следующим образом. В качестве предварительного распределения удельного сопротивления в слое очевидным образом (путем осреднения удельного сопротивления внутри выделенной подобласти) можно задать распределение, показанное на рис. 7, а. В результате одного шага борновской инверсии со сглаживанием было получено распределение удельного электрического сопротивления, показанное на рис. 7, б. Соответствующее ему распределение аномальной проводимости приведено на рис. 7, в.
Из приведенных на рис. 7 результатов видно, что предполагаемый объект находится на том месте, на котором он задан, однако его удельное электрическое сопротивление должно быть ниже, что соответствует истинной геоэлектрической модели.
Рассмотрим теперь ситуацию, когда начальное расположение объекта не совпадает с истинным. Эта ситуация приведена на рис. 8, а. На рис. 8, б приведено распределение удельного электрического сопротивления, полученное в результате одного шага борновской инверсии, а на рис. 8, в - соответствующее ему распределение аномальной проводимости.
-6000 -4000 -2000
2000 4000 6000 8000 10000
-10000 -8000 -6000 4000 -2ООО
2000 4000 6000 8000
Л/с. 7. Восстановление избыточной проводимости при примерном совпадении положений истинного и аппроксимирующего его объектов:
а - предварительное распределение удельного сопротивления в слое; 6 - удельное сопротивление, полученное в результате одного шага борновской инвераш со сглаживанием в выделенных подобластях; в - распределение аномальной проводимости
Из приведенного на рис. 8, в распределения аномальной проводимости с очевидностью следует, что объект в стартовой модели должен быть смещен к левому краю расчетной области (что также соответствует местоположению объекта в истинной модели).
ЗАКЛЮЧЕНИЕ
Разработаны методы решения трехмерных обратных задач геоэлектрики с индукционным источником, основанные на минимизации функционала невязки вдоль направлений, получаемых с помощью борновских приближений, которые могут быть использованы на этапе построения стартовой геоэлектрической модели среды, представляющей собой набор ЗБ-макрообъектов, уточнение параметров которых выполняется в ходе последующего применения процедур нелинейных ЗБ-инверсий, а также для уточнения наличия и местоположения макрообъектов на промежуточных итерациях нелинейной ЗБ-инверсии. С использованием синтетических данных, полученных с помощью ЗБ-моделирования нестационарных геоэлектромагнитных полей для геоэлектрических моделей различной сложности, показана их работоспособность .
-10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 10000 -10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000
а б
-10000 -8000 -6000 -4000 -2000
2000 4000 6000 8000 10000"
О .005 О .004 О .003 О .002 0 .001 0
-0 .001 -0 .002 -0 003 -0 004 -0 005 -0 .006 -0 .007 -0 .008 -0 .009 -0 01 -0 .011 -0 012 "-0 .013
Рис. 8. Восстановление избыточной проводимости при несовпадении положений истинного и аппроксимирующего его объектов:
а - предварительное распределение удельного сопротивления в слое; б - удельное сопротивление, полученное в результате одного шага борновской инверсии со сглаживанием в выделенных подобластях; в - распределение аномальной проводимости
Показано, что применение на заключительных этапах инверсии разработанного метода регуляризации обратной трехмерной задачи, базирующегося на сглаживании удельной проводимости по отдельным подобластям, существенно повышает качество восстановления геоэлектрической модели.
В целом же, разработанные алгоритмы ЗБ-инверсий позволяют в технологиях зондирований с источником в виде токовой петли достаточно уверенно выделять приповерхностные объекты в горизонтально-слоистой среде и определять места дефицита и избытка проводимости относительно среды с уже выделенными ЗБ-объектами.
СПИСОК ЛИТЕРАТУРЫ
[1] Жданов М.С. Теория обратных задач и регуляризации в геофизике / М.С.Жданов. - М.: Научный мир, 2007.-712 с.
[2] Siripunvaraporn W. Three-dimensional magnetotelluric inversion: Data-space method / W. Siripunvaraporn, G. Egbert, Y. Lenbury, M. Uyeshima // Physics of the Earth and Planetary Interiors. - 2005. - 150.- Pp. 3-14.
[3] Gribenko A. Rigorous 3D inversion of marine CSEM data based on the integral equation method / A. Gribenko, M. Zhdanov // Geophysics, 2007. - 72. - № 2. - WA73-WA84.
[4] Cox L.H. 3D inversion of airborne electromagnetic data / L.H. Cox, G.A. Wilson, M.S. Zhdanov // Geophysics. -2012. - 77. - № 4. - P. WB59-WB69.
[5] Персова М.Г. Компьютерное моделирование геоэлектромагнитных полей в трехмерных средах методом конечных элементов / М.Г. Персова, Ю.Г. Соловейчик, Г.М. Тригубович // Физика Земли. - 2011. - № 2. - С. 3-14.
[6] Соловейчик Ю.Г. Метод конечных элементов для решения скалярных и векторных задач / Ю.Г. Соловейчик, М.Э. Рояк, М.Г. Персова - Новосибирск: НГТУ, 2007. - 896 с.
[7] Соловейчик Ю.Г. Моделирование нестационарных электромагнитных полей в трехмерных средах методом конечных элементов / Ю.Г. Соловейчик, М.Э. Рояк, B.C. Моисеев, Г.М. Тригубович // Физика Земли, 1998. - № 10. -С. 78-84.
[8] Badea Е.А. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials / E.A. Badea, M.E. Everett, G.A. Newman, O. Biro // Geophysics, 2001. - Vol. 66. - № 3. - Pp. 786-799.
[9] Могилатов B.C. Математическое обеспечение электроразведки ЗСБ. Система «Подбор» / B.C. Могилатов, А.К. Захаркин, А.В. Злобинский // Новосибирск, АИ «ГЕО». - 2007. - 157 с.
[10] Персова М.Г. О подходе к выполнению ЗО-инверсий данных зондирований нестационарным электромагнитным полем на основе конечноэлементного моделирования / М.Г. Персова, Ю.Г. Соловейчик, Д.В. Вагин, М.Г. Токарева, А.В. Чернышев // Научный вестник НГТУ. - 2011. - № 2 (43). - С. 97-106.
Персова Марина Геннадьевна, доктор технических наук, профессор, профессор кафедры прикладной математики НГТУ. Основное направление научных исследований: конечноэлементное моделирование электромагнитных полей в задачах геоэлектрики и электромеханики. Имеет более 100 публикаций, в том числе 2 монографии. E-mail: persova@fpm.ami.nstu.ru
Симон Евгения Игоревна, м.н.с. кафедры прикладной математики НГТУ. Основное направление научных исследований: решение обратных задач геоэлектроразведки. Имеет 12 публикаций. E-mail: simona_@mail.ru
Соловейчик Юрий Григорьевич, доктор технических наук, профессор, заведующий кафедрой прикладной математики НГТУ. Основное направление научных исследований: конечноэлементное моделирование электромагнитных и тепловых полей. Имеет более 130 публикаций, в том числе 2 монографии. E-mail: kpmt@fpm.ami.nstu.ru
Кошкина Юлия Игоревна, магистрант кафедры прикладной математики НГТУ. Основное направление научных исследований: решение обратных задач геоэлектроразведки. Имеет 2 публикации. E-mail: koshkina_yui@mail.ru
Persova M.G., Simon E.I., Soloveychik Yu.G., Koshkina Yu.I.
Algorithmsfor 3D- inversion ofTEM SoundingData by UsingBornApproximations
The algorithms for 3D- inversion of TEM sounding data based on the inverse linear modelling by Born approximations are considered. The mathematical apparatus and the results of executing these algorithms when using the synthetic analogues of exploration data obtained by mathematical 3D modelling are presented.
Key words: 3D-inversions, Born approximations, 3D-modeling.
