Метод решения обратной задачи магнитотеллурического зондирования по синтезированному магнитному полю Текст научной статьи по специальности «Физика»

УДК 517.958

П.С. Белкин, В.И. Дмитриев

МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОТЕЛЛУРИЧЕСКОГО ЗОНДИРОВАНИЯ ПО СИНТЕЗИРОВАННОМУ МАГНИТНОМУ ПОЛЮ

(кафедра математической физики факультета ВМиК, e-mail: pbelkin@mail.ru)

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

Введение. Метод магнитотеллурического зондирования (МТЗ) основан на использовании измеренного на земной поверхности естественного переменного электромагнитного поля Земли (магнитотеллурического поля) для изучения строения земной коры. Источники МТ-поля находятся на большой высоте в ионосфере и магнитосфере Земли, что приводит к медленности изменения МТ-поля вдоль земной поверхности. Данная особенность МТ-поля позволяет использовать в задачах МТЗ модель плоской волны, нормально падающей на земную поверхность, где электромагнитное поле удовлетворяет уравнениям Гельмгольца:

АЕх (у, z) + шцо(г)Ех (у, z) = 0, АНу (у, z) + шцо(г)Ну (у, z) = 0.

Компоненты электромагнитного поля зависят от характеристик источника (амплитуды плоской волны), которые меняются со временем и не могут быть непосредственно использованы для определения электропроводности земной коры. Для преодоления этого препятствия используются дополнительные краевые условия — соотношения между компонентами электромагнитного поля на земной поверхности. Данные соотношения могут быть записаны в виде

Ех = ZXXHX + ZXy Ну, Еу = ZyxHx + ZyyHy, (1)

где Zxx, Zxy, Zyx, Zyy — компоненты тензора импеданса. Компоненты тензора импеданса зависят только от распределения а (М) электропроводности в Земле и частоты шине зависят от источника поля. При проведении полевых исследований методом МТЗ на земной поверхности многократно измеряются значения тангенциальных компонент МТ-поля, по которым затем находятся компоненты тензора импеданса для набора частот. Для решения обратной задачи (определения электропроводности в Земле) требуется многократно решать прямую задачу (по распределению электропроводности в Земле рассчитывать импеданс на поверхности), постепенно уменьшая невязку между моделируемым импедансом и экспериментальным импедансом. При этом, решая прямую задачу, для расчета импеданса нам необходимо вычислять компоненты как электрического поля, так и магнитного. Это требует большого количества вычислений на каждом шаге итерационного процесса решения обратной задачи. Количество необходимых вычислений можно существенно уменьшить, если изначально восстановить по измеренному импедансу магнитное поле и решать обратную задачу не по импедансу, а по восстановленному магнитному полю.

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

1. Синтез магнитного поля. Рассмотрим далее метод трансформации тензора импеданса в синтезированное магнитное поле. Для этого перепишем линейные соотношения полей на земной поверхности (1) в виде

НХ — YXXEX — YXyEy = 0, Ну — YyXEX — YyyEy = 0, (2)

где

т =

У У

1 XX 1 ху

У У

1ух 1уу:

= [я-1]

тензор адмитанса, вычисляемый через тензор импеданса следующим образом:

2

V _ Яуу

I Т.Т. -

хх — гу2 ' ху эф

V _ У-щ

I Т.11 -

'/■1 ■ ±ух

эф

V _ /х

I ИТ. -

_ у —

'/■1 ' ' ИИ ,/■! '* лгЬ '*

эф

эф

где И2ф = УХХУ„„ — У..,.,,У.,,.,.. Аномальное электрическое поле на земной поверхности можно выразить через магнитное поле

У'

-гш/л о

2тг

— сю сю

'у У-!

гшцо

2тг

Н$°(хо,Уо) йхрйур хо)2 + (у - уо)2 + г2 '

Щ°(хо<,Уо)<1хо<1уо хо)2 + (у - уо)2 + г2

(3)

Подставив (3) в (2), получим

Нах°(х,у)

Щ°(х,у)

— сю сю

7ху(х, у)Н£°(хо, у о) - -ухх(х,у)Ща(хр,уо) ж0)2 + (у - Уо)2

Ъу у)Щ°(хо-, Уо) - ъЛх, у)Ща(ха, у0)

^(х - хо)2 + (у - уо)2

йхойуо = Рх(х,у),

йхойуо = Ру(х,у),

где

7ар(х,у) = г-^-Уа/з(х,у),

Ех(х,у) = ¥хх(х,у)Е? + Уху(х,у)Е* - II).

Ру{х,у) = ГуХ{х,у)Е% + ¥уу(х,у)Е^ - //*.

Полученные интегральные уравнения позволяют определить аномальное магнитное поле и Ну0 на земной поверхности при известных тензоре адмитанса [У] и нормальных полях. В двумерном случае Е-поляризованного поля (Ех(у, г), Ну(у, г), Нг(у, г)) формулы упрощаются и мы имеем

Щ(У,*) =

7Г

Е${у,г) =

П

1 [ Щ°(Уо)(Уо ~ У)йу0

('У ~ Уо)2 + ¿2

Я«°(уо)1п

7Г

— сю сю

\/(|/ - Уо)2 + г2

#0,

Щ(у,г) =

ЖШ^о

Е?(у0) {У° У)2 "\dyo-((у ^ Уо)2 + ^2)

Интегралы понимаются в смысле главного значения.

2. Решение обратной задачи. Рассмотрим решение двумерной обратной задачи с использованием синтезированного магнитного поля. Для решения прямой задачи мы воспользуемся методом интегральных уравнений [2, 3].

В двумерном случае электромагнитное поле имеет две поляризации — ^-поляризацию (Е = = {£^,0,0}, II = {0,1^,0}) и Деполяризацию (Ё = {0, £^,0}, II = {Нх, 0,0}). Для нас представляет

интерес лишь ¿'-поляризация, так как в случае Деполяризации магнитное поле постоянно на земной поверхности и равно удвоенному первичному (падающему) полю.

Рассмотрим двумерную горизонтально-неоднородную модель, описываемую следующим распределением электропроводности:

= \ax(z) при M(y,z)</T, \aT(y,z) при M (у, z) G Т,

где сдг представляет собой электропроводность слоисто-однородной среды, заданную кусочно-постоянной функцией

{со(г) при г < О,

ап при < г < п = 1,..., Ж, оо при г > гдг,

где ^о = 0, — = Кп (кп — толщина п-го слоя); область Т — некоторая конечная двумерная область электропроводности от (у, г) •

Единственную компоненту Ех(у,г) электрического поля представим в виде суммы нормального электрического поля Ех (г), порождаемого сторонними токами в слоистой среде, и аномального электрического поля Ех(у,г), связанного с наличием в слоистой среде проводящей неоднородности Т. Полное поле Ех (у, г) удовлетворяет неоднородному уравнению Гельмгольца

АЕх(у, z) + шцсг(г)Ех(у, z) = -гш/xj.

ст X 1

где А = -щр- + — двумерный оператор Лапласа, ^ — сторонние токи. Нормальное поле удовлетворяет уравнению

•ст IX '

AEx(y,z) + ш[mjjv (z)Ех (у, z) =

Вычитая уравнение (5) из уравнения (4), мы найдем уравнение для аномального поля:

АЕх (у, z) + iwnoN{z)Ex (у, z) = -iwnjx (у, z),

где jx (у, z) рассматривается как избыточный ток в неоднородности Т:

£(У<, z) = (бгг(у, z) - aN (z))Ex (y, z) .

Поля Ex(y,z) и Ex (z) и их первые производные всюду непрерывны, поэтому для аномального поля выполняются условия сопряжения на границах раздела сред, а именно:

(4)

(5)

(6) (7)

[Ех ] =

дЕх дп

= 0,

где п — нормаль к границе разрыва функции сг(у, z). Наличие бесконечной проводимости при г > гдг приводит к граничному условию

= 0.

Z = ZN

Кроме того, функция Ех удовлетворяет условию излучения на бесконечности: Ех убывает как 1 /у/г, где г = у/(у - у0)2 + г0)2 ->• оо.

Для получения интегрального уравнения с целью нахождения полного поля Ех(у,г) необходимо знать Ои — функцию Грина слоисто-однородной среды, зависящую от положения двух точек М(у, г) и М0(уо,го) и определяемую как решение следующей задачи:

AGN(y,z,y0,z0) + iwfiaN(z)GN(у, z,y0, z0) = -шц6(у, z,y0, z0), 'ÔGN'

G (y,z,y0,z0)

= 0,

Z = ZN

(8)

[Gn] =

dn

= 0 на границах раздела сред,

Gn убывает как 1 /у/г, где г = у/(у — уо)2 + {z — zq)2 —> оо.

Здесь 6 (у, г, уо, 2о) = 6(у — уа)8(г — го) — двумерная дельта-функция Дирака, определяемая следующим образом:

Ду, *)%, г, уо, го) ¿,3 = |/(Уо' ^ (9)

М ^ ^ 10, Мо(уо,^А

где /(у, г) — частично непрерывная функция в некоторой области /) и /) /) + ;Ь* — область I) вместе со своей границей 5.

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

поля

. Для того чтобы это сделать, мы выбираем некоторую точку М0(уо,го) и рисуем вокруг нее круг Од радиуса Д, достаточно большого, чтобы область (2, содержащая источники поля, попала внутрь области Од, ограниченной кругом Од, т.е. <3 € Од.

Применим теорему Грина к функциям Ех и С?^ в области Од:

(д^ до^ - ^ ds = f (е?^ - о^) л, (ю)

Од Сд

где ^ обозначает производную по направлению нормали к кругу Од.

Подставив выражения для АЕх и Аиз уравнений (5) и (8) в последнее уравнение (10) и принимая во внимание (9), получим

шцо Л а \Г;.г йз - у - & = гь>ЦоЕ%(у, г)

Я Сд

в случае, если М(у, г) € Од.

Функции С?^ и удовлетворяют условию излучения на бесконечности, следовательно, если мы будем увеличивать радиус окружности Од до бесконечности, то интеграл по Од будет равен нулю, и в результате мы имеем

Е* ^ = II ^

а

Используя аналогичный подход и начиная с уравнения (6), можно получить

Ех(г)=11° (У^Усъ^оЬх (У^) (и)

а

Подставляя выражение (7) в последнее уравнение и учитывая, что Ех(у, г) = Е£ (г) + Е£(у, г), получим интегральное уравнение Фредгольма второго рода для нахождения электрического поля в области неоднородности Т:

Ех(уо, га) = Ех (¿о) + гшц у J [ам(г) - стт(у, г)] О (у, г, у0, га)Ех(у, г) йуйг, (12)

г

где М0(уо, го) € Т.

Если точка М0(уо, ¿о) ^ Т, то уравнение (12) превращается в интегральное выражение, позволяющее рассчитывать электрическое поле Ех(уо, го) в любой точке пространства при известном значении Ех(у, г) внутри неоднородности. Компоненты магнитного поля могут быть рассчитаны по формулам, вытекающим из уравнений Максвелла:

.. ___' "К,- .. г дЕх

шрь аг шрь оу

Или, учитывая (11),

Ну(Уо, = Ну (г0)~— I [ам(г) - ат (у, г)] Ех(у, г) йуйг,

ш/л ,! иг о

i [дом

i Г dGN

Hz(y0, zq) = II:v (z0) + — / — [ctn(z) - aT (y, z)} Ex(y, z) dydz.

J oyo

т

Пусть область неоднородности T представляет собой прямоугольник, который мы разобьем на NM равных прямоугольников Тпт толщины hz и длины hy, п € [1,-N], ш G [I. М\. В каждом прямоугольнике Тпт выберем точку Рпт (Упzm) пересечения его диагоналей. При этом будем считать, что Ну{Рпт) = Ну (уп, zrn), Ех(Рпт) = Ех (уп, zm).

Электрическое поле в неоднородности в данном случае находится из следующей системы уравнений:

n м

Ех{'Уп1 zm) = Ex(zm) + (yi,zj,yn,zm) (aN (z:¡) - aT (y¿, zj)) Ex(yu z:¡)hyhz. (13)

i=i j=i

Система уравнений (13) в матричной форме имеет вид

B = BN + Ga (&N - а) Е, (14)

где

G

[MNxMN] _ О —

/ iw[iGN (yi, Z\, у\, Z\) hyhz

\ílüliGn (yi,zi,yN,zM)hyhz (oT{yi,zi)

vr(y2,Z!)

iw[iGN (yN,ZM,yi,zi)hyhz \ iw[iGN (yN, zM,yN, ZM) hyhz J

a

[MNxMN] _

\

0T(yN,Zl)

0

V

Магнитное поле на земной поверхности определяется по формуле пересчета

N М

НУ (У$> °) = НУ (°)"ЕЕ 01 (аК(гз) - °т (у» гз)) еЛУ»

(15)

г= 1 j=l

где G\ (y,z,y°) =

, y^ — точки земной поверхности q € [1,Q], a Ex(yi,Zj) нахо-

z°=О

дится из системы уравнений (14).

Формула пересчета (15) в матричной форме имеет вид

Н = Ня + С^Ё (а - ам)

(16)

где

G, =

'Gl{yuzl,yl)hvhz ... Gi{yN,zi,yl)hyhz ... Gx (yN, zM, y?) hyhz

KG1(y1,z1,y^)hyhz ... G\ (yjv, zi,yg) hyhz ... Gx (yN, zM, y%) hyhz

(Ex{y bzi) \

Ex (У2; Z\ ) 0

E =

Ex(vn,zi)

a =

0

( o (yi^i) \

v(yi-,z2)

° (yi,zM)

{yNi ZM) J

Ex(yN,zM)J

H =

'Hy{ylo)'

Таким образом, мы определили оператор А решения прямой задачи, позволяющий по известным в области Т значениям электропроводности рассчитать электрическое поле в области Т, а затем магнитное поле на земной поверхности, т. е.

Н = А(ег,ш).

Оператор А является нелинейным по сг, так как электрическое поле в области Т определяется из системы уравнений (14), матрица которой зависит от сг.

Найдем выражение производной Фреше оператора А. Для этого рассмотрим два различных распределения электропроводностей в области Т: а^ и а^р+1\ отличающиеся на некоторую небольшую величину да"'1 = - а"'1. Аналогично Ё(р+1) = Ё!/'' + Ш"' и Н^1) = + <Ш!'"\

Перепишем (16) в виде

н(р+1) = + - = Ня + ¿1 (Ё(Р> + + 6*(р> - <тш

= НМ + аЁМ (<тЬ) _ + 61 ¿Ё^ (а^ - «г*! + аЁ^^) + б^Ё^о-^

Н^ + С, Ж"5

сг

N

61Ё {р)8а{р). (17)

Мы сделали приближение, учитывая, что б^Ё^^сг^ = о (||5сг(р^||). Значение получим, решив систему уравнений (см. (14))

(1 + й0 (<х(р) - о-^)) = -йоЁ^бо-^, = (! + й0 (<х(р) - о-^))'1 (-йоЁ^^^! .

(18)

Из (17) и (18) имеем Н^1) -Н<р> = 61

сг

(р)

сг

N

Г + вп а

г(Р)

сг

N

-1

-60Ё(р)) + Ё^

= Р

+о(|£<гИ|) , (19)

где Р (5сг(р)) — производная Фреше оператора А решения прямой задачи.

Рассмотрим далее решение обратной задачи, которая состоит в определении распределения электропроводности в области Т по найденному экспериментально магнитному полю Нэ. В нашем случае это синтезированное магнитное поле, т. е. магнитное поле, восстановленное по импедансу, измеренному на земной поверхности.

Оператор решения обратной задачи можно записать в виде

а = К (Нэ,ш) .

Обычно найти точное решение обратной задачи практически невозможно. Для решения обратной задачи необходимо найти квазирешение, т. е. такое решение а обратной задачи, что невязка между экспериментальным магнитным полем Нэ и моделируемым полем Н = А(а,ш) будет минимальна. Кроме того, решение должно быть устойчивым к малому изменению начальных данных, т. е. необходимо найти

т£ \г] + аП} , (20)

где

г]= ||НЭ - А{а,ш)\\к =

5 Да; 0

я

2 -1 X \Ну(УЯ1Ш*) ~ Й(уд,Ыа)

5=1 я д=1

(21)

2 И М

- функционал невязки, П = ||сг — сг||ь = X) X) {а{Уг^з) ~ ^{Уг^з)) - стабилизирующий функцио-

1=13=1

нал, а — параметр регуляризации, — частоты, для которых измерено экспериментальное поле, 5 € [1, <5].

Мы выбрали специальную норму функционала невязки для того, чтобы уравнять вклад в магнитное поле на разных частотах, так как используемые частоты могут различаться на несколько порядков. Решение обратной задачи существует, если существует квазирешение. Таким образом, основной задачей при решении обратной является задача минимизации функционала невязки. Эта задача может быть решена точно в случае линейного оператора А решения прямой задачи. В общем случае нелинейного оператора А решение может быть найдено итерационно. Выполним подстановку производной Фреше в функционал невязки (21). Введем обозначения:

= С!

-1

а(р) _ ^ + с0 (&М - ан)) ' (-СоМр)) + Мр) ¿а(р) = 5а +

р = Б1! {5а + ,

нэ _ Н(Р) = г

где ц. — некоторые действительные числа. Перепишем функционал невязки (21):

Г] =

нэ _ н(р+1) = НЭ _ Н(р) _ р1 (6а +

ми

к= 1

в Я . в Я д МИ Я Я д МИ

Е Е1/-1 Е ^ + е ^ (^+

8=1 9=1

8=19=1

Ш,

к= 1

5 Я

А ш

8=19=1

МИ

шя

к= 1

МИ

Е Е -тЧ Е ^ +ф Е + •

8=1 9=1 4 к = 1

к= 1

В силу основной леммы вариационного исчисления для нахождения минимума функционала необходимо, чтобы

—Ф

(¿С

= 0,

с=0

МИ , Я Я ( ми ми

Е н Е Е р™ Е + Е - г>т -

к=1 8=1 д=1 ^ к=1 к= 1

5 Я ( МИ МИ . д Я Я д

Е £К Е + Е =4 = Е Е +^п —

Аш,

шя

V = О,

в=1 д=1

к= 1

к= 1

J

в=1 д=1

о;.«

Таким образом, для определения 5а необходимо решить систему линейных алгебраических уравнений ~К5а = Ь, где

к = Е Е +^Ч = 2 £ х> (ед) ^Ч,

в=1 «=1 Я 8=1 9=1 Я

ь = Е Е (*.Г + *;г) ^ч = 2 £ £ к. (*;г)

в=1 9=1 Я «=1 в=1 Я

Заключение. Мы свели процесс решения обратной задачи к итерационному процессу, на каждом шаге которого решается система линейных алгебраических уравнений. На основании изложенного метода решения был разработан алгоритм и составлена программа.

На рис. 1 представлена модель строения среды, где в нижнем проводящем полупространстве содержится неоднородность ступенчатой формы с тремя областями различной электропроводности. Для этой модели было вычислено магнитное поле //., на земной поверхности {х = 0) для 15 частот. Далее

60 -40 -20 0 20 40 60 Профиль, км

Рис. 3

это поле было взято как экспериментально измеренное и по нему была решена обратная задача описанным выше методом. Так как метод является итерационным, то на рис. 2 приведено решение обратной задачи, полученное на первой итерации. Как видно из рисунка, в этом случае мы получаем некоторую сглаженную область неоднородности. При этом область и по форме и по проводимости существенно отличается от начальной модели. На рис. 3 приведено решение, полученное после десяти итераций. Из рисунка ясно видно, что и размеры и электропроводность решения близки к исходной модели. В дальнейшем к полученному решению можно применить метод контрастирования, позволяющий с помощью использования фильтров более четко выделить зоны различной электропроводности.

СПИСОК ЛИТЕРАТУРЫ

1. Бердичевский М.Н., Дмитриев В. И., Мерщикова H.A. Об обратной задаче зондирования с использованием магнитотеллурических и магнитовариационных данных. М.: МАКС Пресс, 2000.

2. Дмитриев В. И., Б арашков А. С., Белкин П. С. Метод интегральных уравнений в обратных задачах электромагнитных зондирований // Прикладная математика и информатика. № 12. М.: Изд. отдел ф-та ВМиК МГУ, 2002. С. 5-11.

3. Дмитриев В. И., Белкин П. С., Мерщикова H.A. Метод интегральных уравнений в моделировании двумерных задач геоэлектрики // Прикладная математика и информатика. № 18. М.: Изд. отдел ф-та ВМиК МГУ, 2004. С. 5-16.

Поступила в редакцию 24.12.07