CC BY
15
УДК 550.38
DOI: 10.33764/2618-981Х-2019-2-3-47-54
ПРОГРАММНО-АЛГОРИТМИЧЕСКОЕ ОБЕСПЕЧЕНИЕ ТЕХНОЛОГИИ ГЕОМАГНИТНОЙ ТОМОГРАФИИ С ИСПОЛЬЗОВАНИЕМ ДАННЫХ БЕСПИЛОТНОЙ МАГНИТОРАЗВЕДКИ
Михаил Андреевич Максимов
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, научный сотрудник, тел. (905)094-33-22, e-mail: MaksimovMA@ipgg.sbras.ru
Ирина Владимировна Суродина
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Лаврентьева, 6, кандидат физико-математических наук, старший научный сотрудник, тел. (913)903-87-41, e-mail: sur@ommfao1.sscc.ru
Вячеслав Николаевич Глинских
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, доктор физико-математических наук, доцент кафедры геологии месторождений нефти и газа, главный научный сотрудник, тел. (913)910-78-10, e-mail: GlinskikhVN@ipgg.sbras.ru
Реализован высокопроизводительный алгоритм инверсии данных для беспилотного магнитометрического комплекса. Приведены результаты тестирования разработанного алгоритма томографической инверсии данных разновысотной магнитной съемки. Показаны результаты расчетов для тестовых моделей.
Ключевые слова: теория обратных задач, электромагнитная теория, линейные системы, численные решения, БПЛА.
SOFTWARE AND ALGORITHMIC SUPPORT OF GEOMAGNETIC TOMOGRAPHY TECHNOLOGY USING DATA FROM UAV MEASUREMENTS
Mikhail A. Maksimov
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik
Koptyug St., Novosibirsk, 630090, Russia, Researcher, phone: (905)094-33-22, e-mail: MaksimovMA@ipgg.sbras.ru
Irina V. Surodina
Institute of the Computational Mathematics and Mathematical Geophysics SB RAS, 6, Prospect Аkademik Lavrentiev St., Novosibirsk, 630090, Russia, Ph. D., Senior Researcher, phone: (913)903-87-41, e-mail: sur@ommfao1.sscc.ru
Vyacheslav N. Glinskikh
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, D. Sc., Associate Professor of the Department of Geology of Oil and Gas Field, Chief Researcher, phone: (913)910-78-10, e-mail: GlinskikhVN@ipgg.sbras.ru
A data inversion high-performance algorithm is proposed for an unmanned magnetometer complex. We give the testing results of the numerical algorithm for tomographic inversion of data of different-height magnetic survey. We show the results of calculations for the test mode.
Key words: inverse theory, electromagnetic theory, linear systems, numerical solutions, UAV.
Магниторазведка - один из наиболее распространенных методов разведочной геофизики. При необходимости выполнения магнитной съемки в удаленных районах со сложным рельефом оптимальной является многоуровневая съемка с использованием беспилотных летательных аппаратов (БПЛА). Это оборудование позволяет проводить быстрые высокоточные измерения магнитного поля на разных высотах, что дает возможность с приемлемой достоверностью определять параметры геологических объектов таких, как рудные тела, кимберлитовые трубки, а также решать задачи археологии [1]
При таких разведочных работах, когда нет возможности полноценной связи с компьютерными центрами, встает вопрос необходимости быстрых расчетов предварительных результатов инверсии данных с целью определения перспективных участков для более детального исследования, выделения основных аномалий магнитной восприимчивости и оценки их уровня. Для этого необходим быстрый алгоритм инверсии магнитометрических данных с БПЛА, работающий на мобильных платформах, учитывающий рельеф и априорные данные.
Вопросами использования различных математических методов для решения интерпретационных задач магниторазведки занимались, в частности, Блох [2] и Жданов [3]. Ими проведены широкие обзоры большинства методов инверсии и интерпретации и обобщены основные принципы и подходы к решению обратных задач в целом.
Большинство линеаризованных задач сводятся к большим системам линейных алгебраических уравнений (СЛАУ), и существуют различные подходы к ускорению их решения. Наиболее актуальным в современных расчетах является подход параллельных вычислений, позволяющий значительно ускорить вычисление решения больших СЛАУ путем одновременного использования нескольких потоков. В частности, такие возможности предоставляют Intel Open MP (для процессорных вычислительных модулей) и NVIDIA CUDA (для графических ускорителей).
Использование параллельных вычислений при решении СЛАУ для геофизических задач широко освещены в работах Страхова, им предложен метод линейных интегральных представлений. И отдельно стоит отметить его вклад в развитие методики нахождения устойчивых приближенных решений СЛАУ в задачах гравиметрии и магнитометрии [4].
Нами разработан и реализован высокопроизводительный алгоритм решения прямой задачи разновысотной магнитометрии, основанный на магнитоста-тическом подходе к решению уравнений Максвелла, который описан, в частно-
сти, в работе [2] с использованием дискретизации основного уравнения конечно-разностным методом на неравномерной сетке с использованием однородной консервативной схемы [5] и последующей его симметризацией. Данный алгоритм позволяет предсказывать поведение магнитного поля во всем моделируемом пространстве. Это обеспечивает возможность для реализации быстрой инверсии данных беспилотной разновысотной съемки в режиме реального времени. В частности, дополнительной мотивацией использовать подход интегральных уравнений является возможность учета остаточной намагниченности и эффектов перемагничения (для сильномагнитных тел).
В теории интерпретации геофизических полей немаловажную роль играет теория некорректно поставленных задач, классические методы решения которых были предложены в работах [6, 7]. Как известно, обратные задачи геофизики и магнитометрии являются некорректно поставленными, поэтому теории, развитые М.А. Лаврентьевым и А.Н. Тихоновым, успешно применялись для интерпретации потенциальных полей. Соответственно, при реализации алгоритма инверсии особое внимание необходимо уделить вопросам единственности и эквивалентности.
В данной работе представлены результаты тестирования разработанного комбинированного алгоритма трехмерной инверсии магнитных данных (геомагнитная томография). За основу взяты два ключевых подхода: минимизация с использованием градиентов функционала невязки [8] и инверсия линейной системы уравнений путем псевдообращения матрицы чувствительности (SVD) [9].
Традиционно задача инверсии является крайне ресурсоемкой и перед ее реализацией необходимо максимально ускорить алгоритм решения прямой задачи, что и позволяет решение большой СЛАУ на GPU. Чем быстрее работает непосредственно прямая задача, тем более качественную и высокопроизводительную инверсию можно получить за конечное время при прочих равных условиях, так как время решения обратной задачи напрямую зависит от числа определяемых модельных параметров и от времени решения прямой задачи.
Множество алгоритмов инверсии базируется на минимизации функционала невязки различными подходами. Примером такого алгоритма является квази-ньютоновский метод, основанный на изучении поведения градиента и функционала невязки и приближенных вычислений его гессиана. При их применении эффективность вычислений существенно повышается, если имеется возможность быстро и достаточно точно вычислить градиент функционала невязки. Такой возможностью обладает, например, метод сопряженного оператора [10].
Пусть d0 - вектор наблюденных данных магниторазведки с БПЛА (аномальное магнитное поле), dt - вектор сгенерированных синтетических (тестовых) данных путем решения прямой задачи для некоторой модели mt (распре-
деление магнитной восприимчивости в среде), w0 - вектор относительной погрешности наблюденных данных. При решении задачи минимизации одним из часто используемых является функционал невязки следующего вида:
F = l/nYJi(doi -dt,i)2/{w0,id0ii)2,
где n - число измерений, doi - компоненты вектора наблюденных данных, dt i - компоненты вектора текущих тестовых данных, wo i - соответствующие измерениям относительные погрешности. Задача минимизации функционала состоит в нахождении такого набора модельных параметров mt, которая позволит максимально приблизить данные dt к наблюденным d0.
Техника квази-ньютоновской минимизации позволяет быстро и с высокой точностью решить обратную задачу с небольшим числом модельных параметров независимо от априорной информации, поскольку позволяет с высокой вероятностью восстановить глобальный минимум функционала невязки. Это обстоятельство можно использовать для построения референтной модели среды. Мы используем подход Бройдена-Флетчера-Голдфарба-Шанно в условиях ограниченной памяти (L-BFGS), хорошо зарекомендовавший себя в таком классе задач при тестировании различных подходов к минимизации [8].
Далее для референтной среды используется псевдообращение матрицы чувствительности измерений к модельным параметрам, основанное на сингулярном (SVD) разложении. Такой подход позволяет не только быстро и качественно разрешить СЛАУ, исключив линейно-зависимые данные и выделив неразрешимые модельные параметры. Однако для его применения необходимо, чтобы стартовая модель несущественно отличалась от искомой, так как при обращении в качестве решения выдаются лишь небольшие приращения по вектору модельных параметров.
Таким образом, мы имеем комбинированный томографический (линеаризованный) алгоритм, который основан на построении первичной модели низкого разрешения на основе квази-ньютоновского метода и последующего ее уточнения через псевдообращение матрицы чувствительности.
Данный алгоритм реализован в среде Fortran90 с использованием библиотек Intel MKL. На рис. 1 представлена блок-схема реализованного подхода. Синим выделены исходные данные и параметры вычислений. Оранжевым - внутренние вычислительные модули программы. Прямая задача используется каждый раз при расчете итерационных данных и расчете невязки и не указана в цикле повторно для уменьшения громоздкости диаграммы.
Для тестирования возможностей разработанного алгоритма были предложены две базовые модели:
1. Кимберлитовая трубка (цилиндр) под реалистичным рельефом: магнитный цилиндр конечной высоты во вмещающей слабомагнитной среде под пологим рельефом амплитудой 3 м.
Диаметр цилиндра 0=60 м, глубина верхней кромки Но=5 м, высота цилиндра Нс=30м; магнитные восприимчивости цилиндра Кс=5*10-4СИ и вмещающих пород К0=5*10-5СИ.
Задача: определить положение границ по горизонтали и уровень магнитной восприимчивости.
2. Цилиндр и параллелепипед с различной магнитной восприимчивостью и формой:
цилиндр с Б=60м, Н0=5м, Нс=30м, Кс=5*10-4СИ, К0=5*10-5СИ, расположенный Х0=100м, У0=100 м
прямоугольный параллелепипед с Н0=5м, Нр=30м и сторонами 30х100м, Кп=5*10-3СИ, К0=5*10-5СИ., расположенный Х0=-100м, У0=-100м.
Измерения
Грздиент невязки
Модель погрешностей
Функционал невязки
Градиентная минимизация
Рельеф
Итерационные данные
Требуемая точность
Априорная информация
Референтная модель
Стартовая модель
Прямая задача
Тестовые данные
Решение
Рис. 1. Блок-схема программы инверсии
данных разновысотной магнитной съемки с учетом рельефа.
Оранжевым выделены внутренние элементы программы, синим - исходные данные и параметризация задачи
Задача: определить положение тел и относительный уровень магнитной восприимчивости между телами.
Для тестовых моделей были рассчитаны модули вектора магнитной индукции (далее сигналы) на трех высотах: 3м, 5м, 10м. Распределение сигналов в поверхностях полетов для цилиндрической модели показаны на рис. 2.
Алгоритм применялся в двух модификациях: с априорной информацией об общем уровне магнитной восприимчивости в некоторых слоях (фиксированные значения) и без него. Инверсия выполнялась в 3 горизонтальных слоях, каждый из которых был разбит на большое количество малых блоков: в частности, для модели с двумя магнитными телами, при общей сетке разбиения по горизонтали от -240 до +240, размер блока инверсии составил 8 м.
Рис. 2. Расчетные данные для модели «цилиндр» на высотах 3 м, 5м, 10м
и карта дневной поверхности
На рис. 3 представлен результат томографической инверсии данных для цилиндрической модели. Показан слой 2, который включает в себя исследуемое тело. Время расчета с размерами блоков инверсии 12х12м составило порядка 30 минут. Поскольку резких перепадов значений магнитной восприимчивости в полученном решении не должно быть (условие регуляризации решения через псевдообращение), исследуемое магнитное тело несколько размазывается на томографической модели и ее реальную границу можно выбирать по положению среднего значения между относительными максимумами и минимумами на градиентных границах. Из рис. 3 видно, что в такой постановке положение и радиус цилиндра четко определяются в горизонтальной плоскости даже при не очень плотном разбиении на блоки размером 12х12м при общей сетке наблюдений 480х480м. Значение магнитной восприимчивости цилиндра и вмещающей среды определяется с высокой точностью.
Рис. 3. Результат томографической инверсии (слой 2) для цилиндрической модели и расчетные данные на высоте 3м. Размер блоков инверсии 12х12м
На рис. 4 - томографическая инверсия для модели с двумя телами. Основной задачей было восстановить сравнительный и общий уровень магнитной восприимчивости и разделить положение тел в пространстве, что и удалось успешно выполнить. Время расчета для томографических блоков 8х8 м составило порядка 60 минут. Четко выделяются границы обоих тел и правильно определены уровни магнитной восприимчивости.
Рис. 4. Результат томографической инверсии (слой 2) для модели двух магнитных тел и расчетные данные на высоте 3м. Размер блоков инверсии 8х8м
Полученные результаты подтверждают, что предложенная комбинированная схема инверсии с учетом рельефа и априорной информации работает для синтетических моделей и открывает возможность перейти к тестам данных реальных измерений. Использование разновысотных данных позволяет уточнить получаемые результаты за счет, фактически, более высокой точности измерения интерпретируемого поля, изучая характер затухания сигналов, а также корректно учесть влияние рельефа. Открытыми остаются вопросы оптимизации вычислительной схемы для убыстрения расчетов и уменьшения влияния неоднозначности решения по вертикали, возникающеего при применении различных геометрических разбиений пространства томографической инверсии.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Магниторазведчик: нам сверху видно все / М.И. Эпов, А.П. Фирсов, А.В. Савлук, И.Н. Злыгостев // Наука из первых рук. - 2016. - Т. 71-72. - № 5-6. - С. 104-109.
2. Блох Ю.И. Интерпретация гравитационных и магнитных аномалий. Учебное пособие. - М.: РГГРУ, 2009. - 231 с.
3. Жданов М.С. Теория обратных задач и регуляризации в геофизики. - М.: Научный мир, 2007. - 712 с.
4. Страхов В.Н., Страхов А.В. Методы нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенно заданной правой частью, допускающие глубокое распараллеливание вычислений // Вычислительные методы программирования. - 2001, Т. 2, с. 34-55.
5. Самарский А. А. Теория разностных схем. - М.: Наука, 1979. - 655 с.
6. Лаврентьев М.А. О некоторых некорректных задачах математической физики. - Новосибирск: Изд-во СО АН СССР, 1962. - 92 с.
7. Тихонов А.Н. О регуляризации некорректно поставленных задач // ДАН СССР. -1963. - Т. 153. - №1. - С. 49.
8. Fletcher R. Practical methods of optimization. 2nd ed. - New York: John Wiley & Sons, 1987. - 436 p.
9. Стренг Г. Линейная алгебра и ее применения. - М.: Мир.1980. - 456 с.
10. Fichtner A. Full Seismic Waveform Modelling and Inversion. - Berlin and Heidelberg: Springer-Verlag, 2011. - 343 p.
© М. А. Максимов, И. В. Суродина, В. Н. Глинских, 2019
