ФИЗИКА ЗЕМЛИ, 2013, № 3, с. 52-57
УДК 550.837
КВАЗИОДНОМЕРНЫЙ МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ
МАГНИТОТЕЛЛУРИКИ © 2013 г. Н. И. Березина, В. И. Дмитриев, Н. А. Мерщикова
Московский государственный университет им. М.В. Ломоносова, факультет Вычислительной математики и кибернетики, г. Москва E-mail: dmitriev@cs.msu.su Поступила в редакцию 21.05.2012 г.
В статье рассматривается итерационный численный метод решения двумерных обратных задач маг-нитотеллурического зондирования, который дает возможность значительно уменьшить объем вычислительной работы при решении обратных задач в классе квазислоистых моделей. Метод основан на замене оператора прямой двумерной задачи вычисления низкочастотного электромагнитного поля в квазислоистой среде квазиодномерным оператором, который локально аппроксимирует двумерный оператор в каждой точке измерения поля. Метод можно применять при решении обратных задач магнитотеллурики как в случаях E-и ^-поляризаций, так и в случае, когда в обратной задаче одновременно используются значения импедансов для полей обеих поляризаций. Статья содержит описание численного метода и примеры применения рассматриваемого метода к численному решению серии модельных обратных задач магнитотеллурического зондирования.
Ключевые слова: Электромагнитные зондирования, обратные задачи, магнитотеллурика. DOI: 10.7868/S0002333713030046
ВВЕДЕНИЕ
ПОСТАНОВКА ЗАДАЧИ
Работа посвящена квазиодномерному итерационному методу решения обратной задачи маг-нитотеллурического зондирования (МТЗ) в классе двумерных квазислоистых моделей строения среды. Рассматривается обратная задача, которая состоит в определении распределения электропроводности в земной коре по измерениям импеданса магнитотеллурического поля вдоль профиля на земной поверхности в зависимости от координат точки измерения и частоты изменения поля.
Предложенный численный метод решения обратной задачи использует характер распространения электромагнитного поля в проводящих квазислоистых средах. Для квазислоистой модели значение поля в каждой точке наблюдения на профиле в первом приближении можно считать равным значению поля, вычисленному для локальной одномерной слоистой среды, которая получается, если в точке наблюдения провести сечение, перпендикулярное профилю. Это позволяет при решении обратной двумерной задачи значительно уменьшить объем вычислений, за счет того, что в итерационном процессе для большей части итераций численное решение полной прямой двумерной задачи заменяется решением серии прямых одномерных задач.
Рассмотрим двумерную среду, состоящую из заданного числа непересекающихся слоев переменной толщины с постоянными значениями электропроводности в каждом слое и медленно меняющейся в горизонтальном направлении толщиной слоев. Введем систему координат, в которой параметры среды не зависят от координаты x, ось Oz направлена вертикально вниз, в верхнем полупространстве (z < 0) проводимость а = а0 « 0.
Электропроводность а(у, z) в нижнем полупространстве (z < 0) распределена следующим образом:
а (y, z) при z е [0, H], y е [0, l] а = (z) при z е [0,H], y е [0,l] аH = const при z е [H, да).
(1)
Источником электромагнитного поля является плоская волна, зависимость от времени имеет
-Ш г-
вид е , ю — частота колебаний электромагнитного поля.
В двумерном случае система уравнений Максвелла разделяется на две независимые системы: для ^-поляризованного электромагнитного поля с компонентами Е = (Ех,0,0), Н = = (0, Ну, Н .) и ^-поляризованного поля с компонентами Е = (0, Еу, Е.), Н = (Нх,0,0).
В случае Е-поляризации компоненты Ну и Иг магнитного поля выражаются через электрическое поле
1 дих. „ 1 зех
Н =
Н = -— ^, (2)
гюц ду
/юц dz
ц = const — магнитная проницаемость.
Выражение для импеданса электромагнитного поля на земной поверхности (z = 0) в случае E-по-ляризации имеет вид
Ze (у, ю) =
Hy (у, z, ю)
(3)
z=0
В случае Н-поляризации компоненты Еу и Ег электрического поля выражаются через магнитное поле
Еу = -
1 дНх
Ez =
1 дНх
, ^ . (4)
/юц дz * /юц ду
Выражение для импеданса электромагнитного поля на земной поверхности (г = 0) в случае Н-по-ляризации имеет вид
Еу (у, ^ ю)
Zh (У, ю) = —
(5)
z=0
нх ( z, ю)
В данной работе при решении обратной задачи в качестве измеренных данных используется кажущаяся проводимость аа(у, ю), которая вычисляется через импеданс электромагнитного поля по формуле
V (y, т)1
Каждой поляризации поля отвечают свои значения кажущейся проводимости. Для одномерной среды импедансы и кажущаяся проводимость для обеих поляризаций совпадают.
Обратная задача состоит в определении а(у, г) по известным значениям кажущейся проводимости.
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ
Обозначим через N число слоев квазислоистой среды (рис. 1), через а = (ст1, ст2,..., ст я} — проводимость слоев — проводимость подстилающего основания). Предполагается, что верхний слой имеет постоянную толщину, которая принимается равной 1. На земной поверхности в точках у1, у2, ..., уь (точки наблюдения) измерены значения кажущейся проводимости для одной или обеих поляризаций для серии частот колебаний электромагнитного поля ю1, ю2, ..., юм.
Толщина каждого слоя квазислоистой среды зависит от координаты у геоэлектрического разреза. Обозначим через hj = {к(, к{,..., вектор
Л У2 У3 ст0 * 0
_I_I_I___
Уь
СТ1 = 1
! У
^ n =
г
Рис. 1.
толщины слоев геоэлектрического разреза в точке наблюдения, координата которой равна у.
При решении обратной задачи искомыми являются значения проводимости слоев о = {сть ст2, ..., ст^- и/или толщины слоев h = h2, ..., hL} в точках наблюдения.
Так как для одномерной среды импедансы для обеих поляризаций совпадают, квазиодномерный метод решения обратной задачи можно использовать как в случае, когда задан импеданс только для одной поляризации поля (Е или Н), так и в случае, когда измерены значения обоих импедансов.
Решением обратной задачи является распределение электропроводности (значения проводимости слоев о = (ст1, ст2,...,ст я} и/или толщины слоев h = {h1, ..., hL} в точках наблюдения), которые минимизируют целевой функционал невязки.
В случае, когда заданы оба импеданса, целевой функционал представляет собой сумму по всем точкам наблюдений отклонений измеренных значений кажущейся проводимости от значений, полученных при решении обратной задачи.
Ь Ме
Р (°'Ь) = XX \ЛЕ (°'Ь, Ур )- /Е^СУу, Ю т) +
у'=1 т=тЕ
Ь Мн 2
+ XX Лн (Л Ур Ют)- /нЫ(Ур Ют)| ,
у'=1 т=тн
(7)
/>изм / \ 1 ( изм ( \\
/е (У j, ют) = lgio (ств)Е (Уу, ют)),
изм изм
/н (У у, Ют) = lgio (н (у, Ют)) ,
(8)
изм/ \ изм/ \
где а в,Е (уу, ®т) и а в>я( уу, ®т) - измеренные в точке наблюдений у. на профиле значения кажущейся проводимости для Е- и Н-поляризованных полей, ЛЕ (о,И,уу, ют) и Лн (о,И,уу, ют) - операторы прямой двумерной задачи для Е- и Н-поляриза-ций соответственно, а
ЛЕ (о, h, Уу, ®m ) = lgio (®, h, Уу, Ют )) ,
(9)
Лн ( h, У^ Ют ) = 1ё10 (а,Н ( h, ур ©т ))
где а0Е(аДуу, ют) и (o,h,y;, ®т), - значения кажущейся проводимости, полученные в результате решения полных прямых двумерных задач.
Диапазоны частот, которые включаются в целевой функционал, могут не совпадать для Е- и Н-поляризаций. Это дает возможность принимать во внимание различие в поведении Е- и Н-по-ляризованных полей в разных диапазонах частот возбуждающего поля. В случае, когда измеряется только один из импедансов, в выражении (7) для целевого функционала остается только одно из слагаемых, отвечающее тому импедансу, для которого есть данные измерений.
Для численного решения задачи минимизации функционала невязки используется итерационный метод, идея которого состоит в том, что в первом приближении в каждой точке наблюдения оператор прямой двумерной задачи можно заменить на оператор одномерной прямой задачи для локального одномерного разреза
Л1В (с, у}, ют) = ^ю (ааАВ (с, у}, ют)), (10)
где стаШ(с, у], ют) — значение кажущейся проводимости для одномерной слоистой среды, состоящей из N слоев с проводимостью о =
= (а!, а 2,...,а н} и толщиной Ц = {к(, Н{,...,Н^-!}.
На первой итерации вместо функционала (7) минимизируется приближенный функционал невязки
(14)
(15)
(16)
F h) = X X \A1D ( hj, yj, )- fETtij, И.) +
j=1 m=mE (11)
L Mh ()
+ XX Aid ( hj, yj, ®m )- fH3M(yj, ®m)| , j=1 m=mH
в котором вместо операторов AE(о, h,yj, юm) и AH (о, h, yj, ©m) прямой двумерной задачи использованы их локальные одномерные приближения A1D (с, hj, yy, ®m). В результате в каждой точке наблюдения будет получена аппроксимация измеренных данных значениями электромагнитных полей для локальных одномерных разрезов. Затем для двумерного квазислоистого разреза, полученного в результате минимизации приближенного функционала (11), решается полная прямая двумерная задача. Вычисленные при этом значения электромагнитного поля будут отличаться от значений, которые минимизировали невязку для приближенного оператора прямой задачи.
На следующих итерациях в измеренные значения, входящие в функционал (11), вносится поправка, основанная на следующих соображениях.
Для решения {о, h} обратной задачи
AE (о, h, yj, ю„) = frM (yj, fflm), (12)
Ah (о, h, yj, ) = fH3M (yj, ю„). (13)
Локальный одномерный оператор Л1в(<з, h, у, rom) этому же разрезу ставит в соответствие функ-
цию, определяемую локально одномерным импедансом, f (y j, &т) = lgio (ааЛВ (с, h, yj, и„)).
Aw h, yj, Ю т ) = f1 (У j, Ю т ) .
Соотношения
A1D (а h, У i, Ют ) = = flT (yi, ®т) + (( (yi, ®т) - Ae (а, h, y¡, Ют)),
Aid (о, h, y¡, юи) = = /нМ (y j, ®m) + (fl (yj, ®m) - AH (о, h, yj, Юи)) ,
которые следуют из (12)—(14), определяют поправку, которая вносится в измеренные данные для следующей итерации.
На всех следующих итерациях вместо минимизации целевого функционала невязки (11) выполняется поиск минимума следующего функционала
L Me
Fp h) = XX AlD (о, hj, yj, Ют)- fE (yj, ®m )| +
j=1 т=тЕ (17)
L Mh ()
+ XX |AlD (о, hj, yj, Ют)- fH (yj, Ют
j=1 т=тн
где
fE (yj, ®т) = fE'" (j Ют) + + (Г1 (j, Ют) - Af1 (O h, yj, Ют)) ,
fH (yj, Ют) = fH3M (Уу, Ют) + + (Г1 (yj, Ют) - Ah^-1 (o, h, yj, Ют)) , p > 1 — номер итерации.
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.