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

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ иыени М.В. ЛОМОНОСОВА

Р Г 5 ОД

На правах рукописи.

- 6 ЯНВ 191-5

УДК 519.642:537.876.23

Фалалеев Сергей Петрович Математические модели электроыагнитных полей в квазислоистых средах.

Специальность 05.13.16 - применение-вычислительной техники, математического моделирования и математических методов в научных исследованиях.

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук.

МОСКВА 1994

- г -

Работа выполнена в Московском государственном университете имени М.В. Ломоносова на факультете вычислительной математики и кибернетики.

Научный руководитель: доктор физико-математических наук,

профессор В.И. Дмитриев.

Официальные оппоненты: доктор физико-математических наук,

профессор А.Б. Самохин, доктор технических наук, профессор М.Н. Бердичевский.

Ведущая организация: Объединенный институт физики Земли РАН.

Защита состоится " " /р<^р<г\ 199-3'г. в /9

часов минут на заседании диссертационного совета

К.057.05.87 в Московском государственном университете имени М.В. Ломоносова по адресу: 119899 Москва, Воробьевы горы, МГУ, 2-й учебный корпус, факультет вычислительной математики и кибернетики, аудитория (- ¿'£~.

С диссертацией можно ознакомиться в библиотеке факультета ВМиК МГУ во 2 учебном корпусе.

Автореферат разослан " '2С?" декабря 1994' г.

Ученый секретарь диссертационного Совета доцент

В.М. Говоров.

ОБЩЕЕ СОДЕРЖАНИЕ РАБОТЫ.

Актуальность теш. Среди множества методов изучения строения земной коры заметную роль играет метод, основанный на использовании естественного электромагнитного поля, измеренного на земной поверхности, и получивший название метода магнитотеллурического (МТ) зодирования. Задачи определения физических параметров среды по характеристикам поля, распространяющегося в исследуемой среде, являются по своей природе обратными. Эффективность методов их решения во многом определяется качеством используемых методов решения прямой задачи для рассматриваемой модели. Широкое использование персональных компьютеров в прикладных исследованиях активизировало интерес к созданию методов моделирования трехмерных электромагнитных полей в неоднородных средах, некритичных к вычислительным ресурсам.

Классической моделью в исследованиях МТ поля является модель слоистой среды. Применение такой модели правомерно при изучении осадочного чехла с медленным изменением электромагнитных характеристик вдоль земной поверхности. Однако наличие различного рода неоднородностей в земной коре приводит к необходимости рассматривать более сложные трехмерные модели сред, включающие в себя модели неоднородностей того или иного вида. Наиболее широко используемая модель среды представляет собой плоско-параллельную слоистую среду с погруженным в нее трехмерным проводящим телом. Такая модель характерна, например, для задач рудной геофизики, изучения блокового строения земной коры. Методы их решения достаточно громоздки и предъявляют высокие требования к вычислительным ресурсам. Среди множества проблем, представляющих практический интерес, существует широкий класс задач, для которых наиболее естественной моделью является проводящая структура, состоящая из однородных слоев переменной мощности, например, задачи изучения глубинного строения земной коры. Модели такого вида называются квазислоистыми средами. Создание метода моделирования трехмерных электромагнитных полей в таких средах, реализуемого с малыми затратами вычислительных ресурсов, позволило бы существенно продвинуться в решении исключительно сложной трехмерной обратной задачи определения глубинного строения

земной коры.

Цель работы - состояла в создании метода математического моделирования трехмерного магнитотеллурического паля в квазислоистой среде и построении нетребовательного к вычислительным ресурсам, пригодного для использования на персональном компьютере численного алгоритма, реализующего метод.

Научная новизна. В диссертации разработан специальный метод решения трехмерной краевой задачи для системы уравнений Максвелла в квазислоистой проводящей среде. Решение представляется в виде трехмерных электромагнитных потенциалов. Для плотности потенциала получена удобная для численного решения система интегральных уравнений первого рода по поверхностям, разделяющим слои. Предложено новое представление первичного поля в квазислоистой среде, обладающее следующими свойствами: во-первых, построенное поле удовлетворяет, однородной системе уравнений Максвелла во всем пространстве, во-вторых, для построенного первичного поля выполняются условия непрерывности тангенциальных составляющих вне области отклонения границ слоев от нормального разреза и, в-третьих, во всем пространстве вне неоднородности полученное поле совпадает с полем плоской волны, возникающим в слоистой среде с электромагнитными характеристиками нормального разреза. Такое представление позволяет построить экономичный по отношению к памяти компьютера численный алгоритм . решения задачи моделирования электромагнитных полей в квазислоистых средах.

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

Практическая и теоретическая ценность. В диссертации " проведен детальный анализ проблем, возникающих при решении задачи математического—моделирования трехмерного—электромагнитного поля в квазислоистой среде. Предложено ее решение методом интегральных уравнений. Получаемые уравнения максимально учитывают свойства сред и удобны при численной реализации. В ходе решения задачи моделирования предложен оригинальный способ выделения первичного поля, предоставляющий новые алгоритмические возможности. Развитый в диссертации математический метод может применяться при решении задач распространения электромагнитных волн в квазислоистых проводящих средах, возникающих в геофизике и электродинамике.

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

Апробация pa6oW. Содержание работы докладывалось на научной конференции "Теория и практика магнитотеллурического зондирования" ( Москва, 1994 ), научно-исследовательских семинарах лаборатории матемамической физики факультета ВМиК МГУ им. М.В. Ломоносова ( рук. проф. В.И. Дмитриев ), научно-исследовательском семинаре кафедры математической физики факультета ВМиК МГУ им. М.В. Ломоносова (рук. проф. A.M. Денисов ).

Основные результаты диссертации опубликованы в 3 статьях.

Структура и объем работы. Диссертация состоит из введения, трех глав, заключения, списка литературы, включающего 72 названия, и приложения. Общий объем диссертации составляет 117 страниц.

СОДЕРЖАНИЕ ДИССЕРТАЦИИ.

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

В первой главе строится метод решения задачи математического моделирования трехмерных электромагнитных полей в квазислоистых средах. В §1 первой главы описывается модель квазислоистой проводящей среды, в рамках которой - будет исследоваться поведение магнитотеллурического поля. Модель квазислоистой среды определяется как проводящая структура, состоящая из однородных слоев переменной мощности, т.е. электропроводность а может быть представлена в виде :

о(х,у.й) =

о,

о

при йс!^,

ог при (х,у)<2<Ьг+Гг(х,у), (1)

<7Ж1при

Здесь I, 1=1,...,N. номер слоя;

- глубина залегания нижней границы слоя с номером I; Г^(х,у), при 7=1,...,N-1, - непрерывно дифференцируемые функции, обращающиеся в нуль вне области аномалии:

0={(х,у)| |х|<Ьх ,|У|<Ьу }, (2)

10(1,У)нО,

Другими словами, вне локальной области 0 имеем плоскопараллельную слоистую среду со следующим нормальным расцределением электропроводности :

ст(х,у ,й)=аг , . (3)

Завершается §1 постановкой краевой задачи определения электромагнитного поля в квазислоистой среде.

Б §2 первой главы решение поставленной задачи сводится к решению системы поверхностных интегральных.уравнений первого рода с особенностями в ядрах. Основным моментом при применении метода интегральных уравнений для решения краевых задач является подходящий выбор интегрального представления решения. В основе используемого в диссерацш интегрального представления лежит разбиение электромагнитного поля (Е,Н) на сумму двух

о о

полей: первичного поля (Е,Н), содержащего в себе всю информацию об источнике поля и исследуемой среде при отсутствии аномалии и вторичного ( аномального ) поля (Е5,НС), связанного с наличием аномалии в структуре изучаемой модели:

о

Е(х,у,2)=Е(х,у,г) + Еэ(х,у,2), (4)

Н(х,у,а)=Н(х,у,в) +Н3(х,у,2). (5)

Представить юле в виде (4), (5) можно многими способами. В" работе предложен следу кидай способ.

После детального1 рассмотрения вопроса о построении первичного поля в §2 сделан вывод, что обычный в таких случаях прием выделения первичного поля как поля плоской волны в нормальном разрезе приводит к тому, что первичное поле в таком виде не удовлетворяет однородным уравнениям Максвелла. Следовательно, в модель вводятся избыточные токи, что приводит г. интегральным уравнениям, содержащим интегралы 1к. объему, что усложняет метод.__ В диссертации предлагается нетрадиционный способ построения первичного шля, свободный от названного недостатка. Суть его состоит в том, что вне области аномалии О (2) первичное поле полагается равным полю плоской волны, возникающему в слоистой среде с нормальным распределением электропроводности, а в области аномалии внутри каждого слоя аналитически продолжается из регулярного участка слоя. Прояо.гокентте строитср плвлуготм „образом.-_Векторы_первичного___

электромагнитного '3), Нг(х,у,2) в слое с номером" Г,

зададим формулами:

Е (х,у,а)=е

1-0^1е " ~о1°2е )'

(6)

Н'(х,у,2)=п[с;е 1 1 1 + С^е ' 6 ]. где переменная ъ изменяется в пределах:

Здесь I - номер слоя, 1=0,...,N. е,' Л - единичные векторы, задающие поляризацию плоской волны. Константы , соответствующие константы для поля плоской волны в слоистой среде с нормальным распределением электропроводности (3). Для их вычисления известны рекуррентные формулы.

Далее в диссертации предлагается представить Бторичное электромагнитное поле в виде трехмерных электромагнитных потенциалов так, что искомое поле в слое с номер-/,и 7. ;--• 1.....N. имеет вид:

+

Е1(М0)=Ёг(И0) +

Цсг (М.М0)Р-1 (М)ббм+^ёгай|сЗ^Сг(М,М0)7г""1 (!!)№„ + 5г-1 5г-1 ° (?>

Цсг(М,И0)£1(М№м + ^ эг-ай^иц 81(М,М0)£г(М)№м ,

Нг(М0)=Нг(М0) +

л л (8) + 1 |гогм Сггм,м0)5ь1(м)с2зм +1 |г<лм с*гм,м0нг(м№м .

Здесь М - точка с координатами (х,у,г),

М0~ точка с координатами (х0,у0,г0) ,

и - циклическая частота рассматриваемого поля,

р. - магнитная проницаемость исследуемой среда,

0 7 0 7

Е (Мд), Н (М0) - вектор-функции, представляющие собой первичное электромагнитное поле в слое с номером I, введенное ранее,

, - соответственно верхняя и нижняя границы слоя с номером I:

БЬ1={ (х.у.г) 1и=Ьг_1+Г1_1 (х,у)},

Зг={(х,у,2)1г=Ь1+1г(х,у)}.

Поверхностные интегралы в формулах (7), (8) определяют аномальное электрическое и, соответственно, магнитное поле в 1-ом слое, возбуждаемые эффективными токами, распределенными с

плотностью на поверхности и с плотностью {} на

поверхности .

С1, - электромагнитный тензор Грина, который при 1=2,...,N-1 является тензором Грина для однородного пространства с электропроводностью ст^. Он представим в виде:

'10 0 0 1 О О 0 1

С = ц

1к7г е 1

В случае 1=1, С1- тензор Грина пространства с распределением электропроводности вида:

о(гН

-0О , при г<!10, о1 , при г>Ь0.

В случае С" - тензор Грина пространства с распределением электропроводности вида :

гаы , при

Другими словами, первый слой рассматривается вместе с верхним полупространством, а последний - с нижним. Такой выбор тензоров объясняется тем, что поверхности, разделяющие названные области исследуемой среды не содержат аномалий.

Для компактности записи предложенного представления аномального поля в диссертации введены интегральные операторы

Аг(5,5)=1и|ог(М,М0)5(М)с2зм+гГЬ £той[с11им С1 (М.М^еШ) ^

X о

Б Б

|гоГм с1см.и0)е(м)(йы .

(9)

Из условия непрерывности тангенциальных составляющих электромагнитного поля и представления поля (7), (8) в

диссертации получено, что векторы плотности токов ^

удовлетворяют системе интегральных уравнений, которая в операторной записи (9) имеет вид:

- (а1*1^.?) + Аг;1(81+1.ег+1)]]-п11 .

[[вгсзг_1 ) + вЧэ^.е1)] -

- (Вг+1(Б1,?) + В1+1(5г + 1,£1+1))].п1 пг.[[2*(£1(М0)+11(М0)))»пг] =пг.[(н1+1-Н1)-пг] ,

При этом неизвестные векторы лежат в касательной к

поверхности плоскости, поэтому они удовлетворяют уращещям:

(11-пг)=0 , <ег-пг)=0 . •' (12)

Здесь п^ - нормаль к поверхности 1=1.....N-1.

В диссертации отмечены следующие свойства полученной системы:

1) Система (Iи), (И ; является системой интегральны: уравнений первого рода, т.к. матрица свободных членов вырождена. Для ее решения следует применять регуляризирукщий алгоритм.

2) Интегралы, входящие в выражения (9), берутся по

неограниченным поверхностям, то есть векторы , ^ подлежат определению в бесконечных областях. ~

3) Специальный выбор представления первичного поля в (У),

ограничивающих ~ слои7 и" правой частью получаемой системы интегральных уравнений ( правая часть системы отлична от нуля только там, где строение среды- отличается от слоистого ).

4) Равенство нулю правой части интегральных уравнений (10) (11) вне области аномалии 0 вместе с достаточно быстрым убыванием ядер уравнений позволяет при численном решении рассматривать сравнительно малую область интегрирования.

Вторая глава диссертации посвящена построению численного алгоритма решения полученной в первой главе системы поверхностных интегральных уравнений первого рода. В основе разрабатываемого метода лежит метод коллокации с саморегуляризацией за счет выделения особенностей ядер уравнений. В §3 главы II проводится численное исследование вопроса о выборе размеров области интегрирования для аппроксимации интегралов по неограниченным поверхностям. Система интегральных уравнений редуцируется в соответствии с методом коллокации к системе линейных алгебраических уравнений вида:

= К^ДЙЗДЛ-Ь.]]' I ТНЛ «« .^/(^.Т!) ■

(14)

а®сь : 1..1.....Н., .).-!.....Ну, -Г'.

В §4 этой главы рассматривается задача вычисления коэффициентов полученной системы алгебраических уравнений (13), (14). Решаемая система интегральных уравнений (10), (11) является системой первого рода, и для достижения эффекта саморегуляризации требуется находить значения диагональных коэффициентов в алгебраической системе с большой точностью. Существенным моментом при нахождении значений этих коэффициентов системы является то, что при мп=м{ . подынтегральные

функции в выражениях для коэффициентов (10) имеют особенность. Значения требуемых интегралов получаются аналитическим интегрированием особенности. В диссертации подробно рассматривается вычисление всех требуемых интегралов:

1)

7

Подынтегральная функция в интеграле / а (М,М0)с1з имеет

гЛ

«ад

особенность вида 1/Ещ при М-М0 , которая легко интегрируется

о

аналитически;

2) для вычисления интеграла J rot G (M,M)ds используется

------- . Л -----------------------

7

явный вид тензора G , из которого получено, что

Л р t

rot G = rotl-^— + rot g\

где g^- аналитическая тензор-функция. Особенность подынтегральной функции полностью определяется первым слагаемым. В

гк^г

диссертации доказано, что интеграл от rotl-^j- равен нулю,

как интеграл от нечетной функции по симметричной области. Интеграл понимается в смысле главного значения по Коши.

3) вычисление значения выражения grad J dtvм G(MQ,M)dsM -

Q- •

нетривиальная математическая задача, т.к. в него входят производные от интеграла в смысле главного значения по параметру MQ. Для его вычисления в диссертации предлагается оршчшальный-злгорйТм,._состоящий в следующем. Исходя из явного

вида тензора G :

ik7r-1 eiKlv div G = --p-(rx,ry,rz) + div g (p,z+zQ) ,

здесь p= ^/(х-х0)г+(у-у0)2 ,

в диссертации предлагается вычислять значение интересующего выражения как сумму двух слагаемых, а величины каждого из них находить отдельно. Вместо производных по х, у, z от первого слагаемого вычисляются производные по направлениям 1^(1,0,z ), 12=(0,1 ,zy),l3=(zy,zy,~1) . Зная значения этих производных, компоненты градиента легко определить, решив соответствующую систему линейных алгебраических уравнений.

Производная по направлению 11 вычисляется как предел разностного отношения: ¿>F

¿17

|r(x0+5-x,y0-y,z0+zx6-zij+zx(xi-x)+zy(y;j-y) )dsdxdy-

Ql. Ij

■Jf (x0-6-х,y0-y, z0-zxe-ziD.+zx (xi-x)+zy (у .-У) )dsdxdy] -1—].

Qid

{к7г-1 е1Чт Здесь Г=—±--?-(гх,г ,гв) , Р=/ ГОя.

о!. 13

После замены переменной интегрирования х=х'+5 в первом интеграле и х=х'-е во втором интеграле, и вычисления интеграла по х' по формуле среднего, в пределе получено, что уо

_<?р а

у уз-1

- /Г (Х0-Х1, У0-У, (Х^)(у0-у) )йЗ(1у) .

У3'-1 1

В диссертации отмечается, что все равенства точные. Аналогично

находится значение производной по направлению 12.

Производная по направлению 13 аппроксимируется центральной

разностной производной :

¿р^Г /Г(х0-гх5-х,у-2 б-у,20н5-213+йх(х0-х)+2у(у0-у))(1з(1х(Зу -3 7

о'.

13

-¡1 (х0+2хе-х,у+губ-у,204б-21д.+2х(х0-х)+гу(у0-у) )(1Б(1Х(1у] .¿ь\,{

О1.. 13

Таким образом, в диссертации доказано, что интегралы, входящие в выражения для производных, могут быть преобразованы к интегралам, не содержащим особенностей. Численное вычисление последних проводится стандартными средствами. Аналогично вычисляется значение второго слагаемого в рассматриваемом выражении. Отличие лишь в том, что в'качестве направлений 1=1,2,3, используются следующие направления:

11 = (1,0,-2Х),12=(0,1 ,-гу),13=(-ах,-ау.-1) .

В §5 второй главы решается проблема ускорения сходимости и повышения точности расчета значений используемых интегралов Фурье-Бесселя, путем преобразования к сумме двух интегралов, один из которых медленно сходится, но может быть взят аналитически, а другой достаточно быстро сходится и может быть вычислен с помощью стандартных программ.

В §6 второй главы анализируется структура матрицы полученной системы, в которой выделяются отличные от нулевых

блоки. Указывается способ перенумерации, при котором все ненулевые блоки локализуются вблизи главной диагонали. На основе анализа структуры матрицы предлагается блочный итерационный метод решения алгебраической системы (13), (14). Проводится сравнение скорости сходимости' предлагаемого метода с методом простой итерации с оптимальным параметром, блочными методами Якоби и Зейделя.

В §7 второй главы рассматривается алгоритм расчета магнитотеллурических полей на земной поверхности и обсуждается способ применения развиваемого в диссертации метода для расчета используемых в геофизике характеристик среды: тензоры импеданса и кажущегося сопротивления.

Третья глава целиком посвящена вопросам, возникающим при практическом использовании разработанного в первых двух главах подхода к моделированию электромагнитных полей в трехмерных квазислоистых средах. В §8 главы III предлагаются методы повышения эффективности использования ресурсов компьютера. Разработан алгоритм предварительной табуляции значений используемых интегралов Фурье-Бесселя, повышающий скорость вычисления примерно в пять раз при незначительной потере точности. Проводится численное оценивание эффективности применения предлагаемого алгоритма. Рассматривается прием, который при наличии двух плоскостей симметрии у конфигураций аномалий в строении среды позволяет сократить размер матрицы решаемой системы в шестнадцать раз и почти в четыре раза повысить скорость вычисления.

В §9 главы III описывается структура пакета программ, реализующего разработанный алгоритм для фундаментальной модели МТ зондирования, и приводятся некоторые его характеристики. В §10 проводится численное исследование разрешающей способности метода в задаче геометрического зондирования квазислоистых сред. Моделирование частотного зондирования с помощью созданного программного комплекса проводится в §11 третьей главы диссерации. Исследуются асимптотические свойства магнитотеллурических полей при неограниченном росте протяженности горизонтальных размеров аномалии вдоль одной из осей.

В заключении формулируются основные результаты, полученные в диссертации, состоящие в следующем:

1) Разработан специальный метод решения трехмерной краевой задачи для системы уравнений Максвелла в квазислоистой проводящей среде. Решение представляется в виде трехмерных электромагнитных потенциалов. Для плотности потенциала получена удобная для численного решения система интегральных уравнений первого рода, где интегралы берутся по поверхностям, разделяющим слои;

2) Построен численный метод решения полученной системы интегральных уравнений первого рода, использующий эффект саморегуляризации и не предъявляющий жестких требований к вычислительным ресурсам. Алгоритм адаптирован для реализации на персональном компьютере;

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

Публикации.

Основные результаты диссертации опубликованы в работах:

1. Дмитриев В.И., Фалалеев С.П. Метод трехмерных потенциалов в задаче моделирования электромагнитных полей в квазислоистых средах. // Вестник МГУ, сер. 15, Вычислительная математика и кибернетика. 1993 г., М, с.30-40.

2. Дмитриев В.И., Фалалеев С.П. . Численное исследование магнитотеллурического поля в квазислоистой среде. // Изв. АН СССР, Физика Земли. 1993 г., №11, с.55-60.

3. Дмитриев В.И., Фалалеев С.П. Математическое моделирование электромагнитных полей в квазислоистых средах. // Математическое моделирование и решение обратных задач математической физики. М. Изд. МГУ 1994 г., с.170-178