S e MR ISSN 1813-3304

СИБИРСКИЕ ЭЛЕКТРОННЫЕ

МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ

Siberian Electronic Mathematical Reports

http://semr.math.nsc.ru

Том 5, стр. 620–631 (2008) УДК 517.988.68

MSC 65R30

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ

МАГНИТОМЕТРИИ

В. В. ВАСИН, Е. Н. АКИМОВА, Г. Я. ПЕРЕСТОРОНИНА, П. С. МАРТЫШКО, В. А. ПЬЯНКОВ Abstract. The three-dimensional inverse magnetometry problem is investigated. The magnetometry equation is reduced to the two-dimensional nonlinear equations of the rst kind. For solving nonlinear magnetometry equation the iterative regularized Newton method and the modied Levenberg-Marquardt method are used. Results of the solution of magnetometry problem for real magnetic elds are presented.

Keywords: inverse magnitometry problem, integral equation, parallel regular algorithms, parallel computing system MVS-1000.

1. Введение Основное направление исследований, которое развивалось в рамках интеграционного проекта в Институте математики и механики УрО РАН совместно с Институтом геофизики УрО РАН, было связано с созданием вычислительных технологий решения обратных задач гравиметрии и магнитометрии. Эта технология включает все необходимые этапы: 1) полевые измерения магнитного поля; 2) предварительная обработка измеренного поля и выделение аномального эффекта от поверхности раздела сред с постоянной вертикальной компонентой вектора намагниченности; 3) решение нелинейного Vasin, V.V., Akimova, E.N., Perestoronina, G.Ya., Martyshko, P.S., Pyankov, V.A., Methods for inverse magnitometry problem solving.

c 2008 Васин В.В, Акимова Е.Н., Пересторонина Г.Я, Мартышко П.С., Пьянков В.А.

Итоговый научный отчет по междисциплинарному интеграционному проекту СО РАН:

Разработка теории и вычислительной технологии решения обратных и экстремальных задач с приложением в математической физике и гравимагниторазведке.

Работа поддержана УрО РАН (Интеграционный проект с СО РАН и ДВО РАН) и частично РФФИ (грант 06-01-00116).

Поступила 17 июля 2008 г., опубликована 27 ноября 2008 г.

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОМЕТРИИ

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

За отчетный период сотрудниками ИММ УрО РАН получены следующие основные результаты:

а) разработаны и программно реализованы новые алгоритмы предварительной обработки геофизических полей;

б) итеративно регуляризованный метод Ньютона и метод Левенберга–Марквардта были адаптированы к решению уравнения магнитометрии;

в) все методы были протестированы на квазимодельных данных и использованы для решения серии задач магнитометрии с реальными полями.

Сотрудниками Института геофизики УрО РАН:

а) проведены площадные измерения магнитного поля для трех районов;

б) выполнена предварительная обработка по выделению аномального поля для полученных данных;

в) детально исследована структура аномалиеобразующего объекта в районе Башкирской магнитной аномалии.

2. Постановка задачи Рассматривается трехмерная структурная обратная задача магнитометрии по численному восстановлению разделяющей поверхности сред (геологической границы) на основе данных о магнитном поле, измеренном на некоторой площади земной поверхности, и скачке вектора намагниченности.

Функция z = z(x, y), описывающая искомую поверхность раздела, удовлетворяет нелинейному двумерному интегральному уравнению Фредгольма первого рода bd z(x, y ) (1) B[z] J [(x x )2 + (y y )2 + z 2 (x, y )]3/ ac H dx dy = G(x, y), [(x x )2 + (y y )2 + H 2 ]3/ где J скачок вертикальной компоненты вектора намагниченности, G(x, y) аномальное магнитное поле, обусловленное отклонением искомой поверхности от асимптотической плоскости z = H.

Предварительная обработка магнитных данных, связанная с выделением аномального поля (т.е. получением правой части уравнения (1), выполнялась по методике, разработанной П.С. Мартышко и И.Л. Пруткиным [1].

3. Методы решения базового уравнения После дискретизации уравнения (1) на сетке n = M N, где задана правая часть G(x, y), и аппроксимации интегрального оператора B по квадратурным формулам имеем систему нелинейных уравнений (2) An [z] = Fn.

Для решения системы уравнений вида (2) используется итеративно регуляризованный метод Ньютона [2] z k+1 = z k [An (z k ) + k I]1 (An (z k ) + k z k Fn ) (3) либо метод Левенберга–Марквардта [3] [4], 622 В.В. ВАСИН, Е.Н. АКИМОВА, Г.Я. ПЕРЕСТОРОНИНА и др.

(4) Здесь An (z k ) и Fn конечномерные аппроксимации интегрального оператора и правой части в уравнении (1), An (z k ) производная Фреше оператора B в точке z k, I единичный оператор, k последовательность положительных параметров регуляризации.

Нахождение очередного приближения z k+1 по найденному z k сводится к решению СЛАУ где Ak = An (z k ) + k I несимметричная заполненная n n матрица, Fn = An z (An (z ) + k z k Fn ) вектор размерности n.

Таким образом, для нахождения очередного приближения z k+1 метода Ньютона (3) необходимо решать СЛАУ (5) с несимметричной полностью заполненной n n матрицей, а для нахождения очередного приближения z k+ метода Левенберга–Марквардта (4) необходимо решать СЛАУ с симметричной n n матрицей.

На каждом шаге метода Левенберга–Марквардта используется метод квадратного корня, на каждом шаге метода Ньютона используются прямые методы Гаусса или Гаусса-Жордана или итерационные методы градиентного типа в регуляризованном варианте для решения СЛАУ с симметричной матрицей.

Для реализации итерационных процессов (3) и (4) при решении задачи магнитометрии (1) применялись (с некоторыми модификациями) численные процедуры, описанные для задачи гравиметрии [5] [6].

Численная реализация и распараллеливание итерационных методов выполнены на многопроцессорном вычислительном комплексе МВС–1000/ с помощью библиотеки MPI [7] на языке Фортран.

4. Численные результаты для реальных данных На многопроцессорном вычислительном комплексе МВС-1000/32 были решены следующие задачи магнитометрии с реальными данными [8] [11].

Предварительная обработка измеренных магнитных полей для исследуемых районов, связанная с выделением аномальных полей (т.е. получение правой части уравнения (1), была выполнена сотрудником Института геофизики УрО РАН В.А. Пьянковым.

Задача 1. Для одного объекта, расположенного в Северном Казахстане, были обработаны данные магнитного поля, измеренные на площади 1950 (м2 ) с шагом x = y = 20 (м). Расстояние до асимптотической плоскости составляло H = 200 (м). Скачок намагниченности принимался равным J = 20 (A/м).

В случае дискретизации с шагом сетки x = y = 40 (м) исходная задача сводится к СЛАУ с несимметричной матрицей порядка 1830 1830. В случае использования сетки с шагом x = 20 (м), y = 40 (м) необходимо решать СЛАУ 3600-го порядка.

Задача магнитометрии решалась итеративно регуляризованным методом Ньютона с числом итераций Nн = 5 и k = 1. На каждом шаге базового

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОМЕТРИИ

процесса использовался параллельный метод сопряженных градиентов с числом итераций Nмсг = 15. В ходе решения задачи относительная норма невязки уменьшилась в 10 раз.

На рис. 1 представлены изолинии аномального магнитного поля.

Рис. 1. Изолинии аномального магнитного поля для поверхности S На рис. 2 (шаги сетки x=y=40) представлена восстановленная поверхность раздела S1.

Рис. 2. Восстановленная поверхность раздела S1 (Казахстан) Задача 2. Для рудного объекта был обработан массив магнитных данных Оренбургской аномалии, измеренный на площади 125 147.4 (км2 ) с шагом x = 1.25, y = 2.2 (км). Расстояние до асимптотической плоскости составляло H = 6 (км). Скачок намагниченности принимался равным J = 2 (A/м).

После дискретизации исходного уравнения на сетке задача сводится к СЛАУ с несимметричной матрицей порядка 6700 6700. Задача решалась итеративно регуляризованным методом Ньютона с числом итераций Nн = 7 и k = 0.5 На каждом шаге метода Ньютона использовались параллельные метод простой итерации (Nмпи = 5000) либо метод сопряженных градиентов (Nмсг = 5). В ходе решения задачи относительная норма невязки уменьшилась в 70 раз.

624 В.В. ВАСИН, Е.Н. АКИМОВА, Г.Я. ПЕРЕСТОРОНИНА и др.

Заметим, что оба итерационных метода (МПИ и МСГ) при подходящем выборе параметров регуляризации k дают близкие результаты, что говорит о хорошем качестве решения.

На рис. 3 показана восстановленная поверхность раздела S2.

Рис. 3. Восстановленная поверхность раздела S2 (Оренбург) Задача 3. Для рудного объекта был обработан массив магнитных данных Башкирского заповедника, измеренного на площади 93.75 165 (км2 ) с шагом x = 1.875, y = 3.3 (км). Расстояние до асимптотической плоскости составляло H = 10 (км). Скачок намагниченности принимался равным J = 1 и J = 0.1 (A/м).

После дискретизации исходного уравнения на сетке задача сводится к СЛАУ с несимметричной матрицей порядка 2500 2500. Задача решалась итеративно регуляризованным методом Ньютона с числом итераций Nн = и параметром регуляризации k = 0.01 и методом Левенберга–Марквардта с числом итераций Nл-м = 10 и параметрами k = 50, = 1. На каждом шаге метода Левенберга–Марквардта использовался метод квадратного корня, а на каждом шаге метода Ньютона использовался параллельный метод сопряженных градиентов с числом итераций Nмсг = 4. В ходе решения задачи с начальным приближением z0 (x, y) = 0.3 относительная норма невязки уменьшилась в 600 раз.

На рис. 4 представлена восстановленная поверхность раздела S3 для случая J = 0.1 (А/м), полученная методом Ньютона.

На рис. 5 представлена восстановленная поверхность раздела S3 для случая J = 0.1 (А/м), полученная методом Левенберга–Марквардта.

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОМЕТРИИ

Рис. 4. Восстановленная поверхность раздела S3 (метод Ньютона) Рис. 5. Восстановленная поверхность раздела S3 (метод Итоговые результаты переданы специалистам по прикладной геофизике для геологической интерпретации.

Вывод. Найденные поверхности раздела согласуются с представлением геофизиков об исследуемых районах Казахстана, Оренбурга и Башкирии.

Для методов МСГ и МПИ проведено сравнение коэффициентов ускорения и эффективности последовательных и параллельных алгоритмов. Рассмотрим коэффициенты ускорения и эффективности где Tm время выполнения параллельного алгоритма на MBC–1000 с числом процессоров m (m > 1), T1 время выполнения последовательного алгоритма на одном процессоре.

Tm представляет собой совокупность чистого времени счета и накладных расходов на межпроцессорные обмены 626 В.В. ВАСИН, Е.Н. АКИМОВА, Г.Я. ПЕРЕСТОРОНИНА и др.

Число процессоров m соответствует упомянутому разбиения векторов на m частей.

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

В табл. 1 приведены времена счета на МВС-1000/32 и коэффициенты ускорения и эффективности решения нелинейной задачи магнитометрии (Казахстан) о восстановлении поверхности раздела итеративно регуляризованным методом Ньютона с использованием параллельного МСГ (матрица 3600 3600).

Таблица 1. Метод Ньютона+Метод сопряженных градиентов В табл. 2 3 приведены времена счета на МВС-1000/32 и коэффициенты ускорения и эффективности решения нелинейной задачи магнитометрии (Оренбург) о восстановлении поверхности раздела итеративно регуляризованным методом Ньютона с использованием параллельных МСГ и МПИ (матрица 6700 6700).

В табл. 4 приведены времена счета на МВС-1000/32 и коэффициенты ускорения и эффективности решения нелинейной задачи магнитометрии (Башкирия) о восстановлении поверхности раздела итеративно регуляризованным методом Ньютона с использованием параллельного МСГ (матрица 2500 2500).

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

5. Исследование структуры аномалиеобразующего объекта в При построении магнитной модели, как и любой другой модели земной коры, немаловажным фактором является то, какая геодинамическая концепция заложена в ее основу. Известен тот факт, что Уральские структуры представляют собой пространственно стабильный более чем 2000-километровый линеамент, представленный рядом массивов (Денежкин Камень, Кумба, Кытлым и т.д.). Эти крупные массивы располагаются в центральной части Урала и пространственно разделены. Поскольку

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОМЕТРИИ

Таблица 2. Метод Ньютона+Метод сопряженных градиентов Таблица 3. Метод Ньютона+Метод простой итерации (6700 6700) Таблица 4. Метод Ньютона+Метод сопряженных градиентов деформирование горных пород в геологическом масштабе времени можно рассматривать как следствие движений сильно вязкой среды, то такое размещение ультрамафитовых массивов находит свое объяснение [12].

Подобные пространственные закономерности характерны для развития диапира. Неотъемлемой особенностью формирования диапировых структур является и то, что на периферии диапира должны возникать более плоские локальные структуры антиклинального типа. Отражение реликтов такого явления наблюдается наиболее ярко в региональном гравитационном поле Урала и проявляется в виде трехгорбых положительных структур поля 628 В.В. ВАСИН, Е.Н. АКИМОВА, Г.Я. ПЕРЕСТОРОНИНА и др.

g. Достаточно контрастно такие особенности развития структур Урала проявляются и в электрической модели земной коры. Анализ комплекса данных позволяет подойти к построению магнитной модели земной коры вполне подготовленными.

Известно, что обратная задача магнитометрии для однородных односвязных областей однопараметрически неоднозначна по параметру I (определяется магнитная масса Mx ). То есть существуют целые однопараметрические семейства эквивалентных областей. Показано, что для обратных задач структурной магнитометрии имеет место трехпараметрическая неоднозначность по параметрам: |I|,, h (I скачок интенсивности намагниченности, направление намагниченности, h положение асимптоты). Таким образом, для выбора вариантов строения среды необходимо получение единственного решения обратной задачи [12].

В связи с возникновением в последние десятилетия высокоточной и стабильной аппаратуры нового класса квантовых магнитометров появилась возможность получения дополнительной информации, необходимой для определения единственного решения обратной задачи магнитометрии. Так, например, комплексная интерпретация временного аналога профильной магнитовариационной кривой, полученной во вращающемся поле Sq – вариаций, и аномальной магнитной кривой позволяет получить единственное решение обратной задачи [13].

Для комплексной интерпретации статического магнитного поля и аномальных вековых вариаций на профиле Шафраново-Магнитогорск (рис. 6) используются данные повторных высокоточных наблюдений на Башкирском полигоне, включающие в себя до десяти циклов измерений за период в несколько лет. На региональном субширотном профиле выделена отрицательная магнитная аномалия с амплитудой в экстремуме до 300 нТл.

По результатам съемки получены изменения разностей T между пунктами векового хода и магнитовариационной станцией Нугуш [14].

Рис. 6. Магнитное поле T и разрез земной коры по профилю Шафраново – Магнитогорск (цифры на линиях в нижней части рисунка –

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОМЕТРИИ

Нами в качестве примера рассматриваются данные аномальных вековых вариаций T по профилю Кага – Старосубхангулово, пространственно совпадающему с отрицательной магнитной аномалии.

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

С целью изучения динамики аномального векового хода и разделения полей от различных источников нами проведен факторный анализ временных рядов T (t). В результате анализа получены две факторные кривые, по своей морфологии напоминающие магнитное поле от объекта с современной намагниченностью и более глубинного объекта, намагниченного против поля, соответственно. Следует учесть, что в отличие от метода аппроксимации сингулярными источниками нет необходимости фиксировать направление намагничивания как для источника положительного, так и отрицательного поля (рис. 7).

Из совместного анализа статического магнитного поля и пространственновременных особенностей вековых вариаций T следует, что аномалиеобразующий объект представляет собой антиклиналь, породы которой намагничены против современного поля. Единственное отличие в том, что региональная кривая векового хода смещена на 5 км относительно соответствующей статической. В результате такой асимметрии удается получить единственное решение для статического поля. Таким образом, региональный объект характеризуется остаточной намагниченность 1 А/м с направлением против современного поля, а его поверхность в экстремальной точке практически касается поверхности земли.

Рис. 7. Разделение аномального векового хода на составляющие от Таким образом, не задавая ни параметров, ни типа источников в случае факторного анализа удалось получить результат аналогичный тому, что получен при автоматизированном подборе методом СИГМА. Более того, в результате анализа временной динамики факторных нагрузок удалось показать, что в случае положительного источника аномального векового хода 630 В.В. ВАСИН, Е.Н. АКИМОВА, Г.Я. ПЕРЕСТОРОНИНА и др.

геомагнитного поля Т его магнитная масса М уменьшается со временем, а магнитная масса М отрицательного источника растет по абсолютной величине.

Таким образом, изучаемый региональный объект действительно намагничен против современного поля, что не является неким исключением для складчатого Урала. Следовательно, комплексирование аномальных вековых вариаций и статического поля при интерпретации магнитных аномалий существенно облегчает решение обратной задачи магнитометрии в рамках трехпараметрической неоднозначности. Имея набор исходных параметров, мы получили решение обратной задачи магнитометрии для трехмерной поверхности (рис. 8).

Рис. 8. Поверхность антиклинальной складки (по осям координат – расстояние и глубина до поверхности объекта в км, соответственно) [1] Мартышко П.С., Пруткин И.Л., Технология разделения источников гравитационного поля по глубине, Геофизический журнал, 25: 3 (2003), 159–168.

[2] Bakushinsky, A., Goncharsky, A., Ill-Posed Problems: Theory and Applications, Cluver Publishers, 1994.

[3] Васин В.В., Мокрушин А.А., Итерационные процессы типа Гаусса-Ньютона для некорректных операторных уравнений, ДАН, 371: 1 (2000).

[4] Васин В.В., Еремин И.И., Операторы и итерационные процессы Фейеровского типа.

Теория и приложения, Екатеринбург: УрО РАН. 2005, 210 с.

[5] Акимова Е.Н., Скорик Г.Г., Регулярные методы решения обратной задачи гравиметрии, Сибирские Электронные Математические Известия, 5 (2008), 509–517.

[6] Акимова Е.Н., Васин В.В., Пересторонина Г.Я., Тимерханова Л.Ю., Мартышко П.С., Кокшаров Д.Е., О регулярных методах решения обратных задач гравиметрии на многопроцессорном вычислительном комплексе, Вычислительные методы и программирование, Москва: МГУ, 8: 1 (2007), 107–116.

[7] Baranov A.V., Latsis A.O., Sazhin C.V., Khramtsov M.Yu., The MVS-1000 System User’s Guide // http://parallel.ru/mvs/user.html.

[8] Akimova E.N., Vasin V.V. Stable parallel algorithms for solving the inverse gravimetry and magnetometry problems, International Journal Engineering Modelling – University of Split, Croatia, 17: 1–2 (2004), 13–19.

[9] Акимова Е.Н., Васин В.В., Пересторонина Г.Я., Решение обратной задачи магнитометрии с использованием параллельных технологий, Материалы 33-й сессии Международного семинара им. Д.Г. Успенского Вопросы теории и практики

МЕТОДЫ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МАГНИТОМЕТРИИ

геологической интерпретации гравитационных, магнитных и электрических полей, Екатеринбург: ИГФ УрО РАН, 2006, 447–450.

[10] Акимова Е.Н., Васин В.В., Решение обратной задачи магнитометрии на многопроцессорном вычислительном комплексе МВС–1000, Материалы 34-й сессии Международного семинара им. Д.Г. Успенского "Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей". – Москва: ИФЗ РАН, 2007, 8–11.

[11] Акимова Е.Н., Васин В.В., Скорик Г.Г., Решение обратных задач магнитометрии и гравиметрии о восстановлении разделяющей поверхности сред, Материалы 35-й сессии международного семинара им. Д.Г. Успенского Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, Ухта: УГТУ, 2008, 8–11.

[12] Федорова Н.В., Цирульский А.В., К вопросу о разрешимости обратной задачи логарифмического потенциала для контактной поверхности в конечном виде, Изв.

АН СССР, Физика Земли, 10 (1976), 61–72.

[13] Мартышко П.С., Пьянков В.А., О единственности решения обратной задачи теории потенциала в методе подмагничивания вращающимся полем Sq – вариаций, Вопросы теории и результаты применения методов интерпретации и моделирования геофизических полей. Свердловск: УрО РАН СССР, 1989, 13–17.

[14] Шапиро В.А., Пьянков В.А., Токовая аномалия векового хода магнитного поля T в Башкирии, Геомагнетизм и аэрономия, 5 (1976), 943–946.

Владимир Васильевич Васин, Елена Николаевна Акимова, Галина Яковлевна Пересторонина Институт математики и механики УрО РАН, ул. С. Ковалевской 16, 620219, Екатеринбург, Россия E-mail address: [email protected], [email protected] Петр Сергеевич Мартышко, Валентин Александрович Пьянков Институт геофизики УрО РАН, ул. Амундсена 100, 620016, Екатеринбург, Россия E-mail address: [email protected], [email protected]