CC BY
10
УДК 550.38
DOI: 10.18303/2618-981 X-2018-3 -241 -247
РАЗРАБОТКА ПРОГРАММНО-АЛГОРИТМИЧЕСКОЙ БАЗЫ ДЛЯ ТЕХНОЛОГИИ ГЕОМАГНИТНОЙ ТОМОГРАФИИ С ИСПОЛЬЗОВАНИЕМ ДАННЫХ БЕСПИЛОТНОЙ РАЗНОВЫСОТНОЙ МАГНИТОРАЗВЕДКИ
Михаил Андреевич Максимов
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 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, доктор физико-математических наук, зав. лабораторией скважинной геофизики; Новосибирский национальный исследовательский государственный университет, 630090, Россия, г. Новосибирск, ул. Пирогова, 2, доцент кафедры геологии месторождений нефти и газа; Новосибирский государственный технический университет, 630073, Россия, г. Новосибирск, пр. К. Маркса, 20, профессор кафедры геофизических систем, тел. (383)330-45-05, e-mail: GlinskikhVN@ipgg.sbras.ru
Выполнена программная реализация алгоритма трехмерной прямой задачи разновы-сотной магнитной съемки с учетом рельефа дневной поверхности с использованием высокопроизводительных вычислений на графическом ускорителе. Реализована принципиальная схема решения обратной задачи с использованием квазиньютоновской минимизации функционала невязки с использованием техники сопряженного оператора для ускорения вычисления градиентов. Проведены тесты на синтетических моделях с небольшим количеством модельных параметров.
Ключевые слова: теория обратных задач, электромагнитная теория, магнитное поле, численные решения, БПЛА.
DEVELOPMENT OF SOFTWARE AND ALGORITHMIC BASES
FOR THE GEOMAGNETIC TOMOGRAPHY TECHNOLOGY USING UAV MAGNETIC
DATA FROM DIFFERENT HEIGHTS
Mikhail A. Maksimov
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, Researcher, Laboratory of Borehole Geophysics, 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 Akademik Koptyug St., Novosibirsk, 630090, Russia, D. Sc., Head of Laboratory of Borehole Geophysics; Novosibirsk National Research State University, 2, Pirogova St., Novosibirsk, 630073, Russia, Associate Professor Department of Geology of Oil-and-Gas Field; Novosibirsk State Technical University, 20, Prospect K. Marx St., Novosibirsk, 630073, Russia, Professor, Department of Geophysical Systems, phone: (383)330-45-05, e-mail: GlinskikhVN@ipgg.sbras.ru
The software implementation of the algorithm for the three-dimensional direct problem of a different altitude magnetic survey is provided taking into account the relief of the day surface with the use of high-performance computations on a graphic accelerator. A schematic diagram of the solution of the inverse problem is realized using the quasi-Newton minimization of the residual functional using the technique of the adjoint operator to accelerate the computation of the gradients. Tests are carried out on synthetic models with a small number of model parameters.
Key words: inverse theory, electromagnetic theory, magnetic field, numerical solutions,
Магнитометрическая разведка, которая дает информацию о локальных изменениях магнитного поля Земли, позволяет сузить круг поисков геологических объектов в труднодоступных регионах страны. Современные методы геомагнитных исследований имеют серьезные недостатки, главные из которых -чрезвычайно высокая стоимость работ и недостаточная точность данных. Актуальной задачей настоящего времени является повышение качества результатов инверсии и интерпретации данных магнитной съемки. Одним из наиболее перспективных подходов является использование разновысотных данных с целью повышения достоверности получаемой информации о геометрии исследуемой среды: в первую очередь анализ глубин залегания различных тел, что при классической наземной или аэросъемке затруднено из-за отсутствия зондирующего параметра на глубину. В ИНГГ СО РАН создан принципиально новый аэромагнитный комплекс на основе беспилотного гексакоптера, оборудованного высокоточным магнитометрическим датчиком. Его применение особенно актуально при картировании больших, практически непроходимых территорий, покрытых сверху водоемами и болотами, а также лесных массивов и зон вечной мерзлоты.
Физико-математические основы для решения прямой трехмерной задачи магниторазведки изложены в работе [1]: используются классические уравнения Максвелла в постановке магнитостатики. Магнитное поле описывается двумя основными характеристиками: напряженностью Н и индукцией В, между которыми существует соотношение:
где |ы() - магнитная постоянная, / - намагниченность в изучаемой точке.
Намагниченность вещества сложным образом зависит от его магнитных свойств и намагничивающего поля, что характеризуется петлей гистерезиса. В магниторазведке применяют упрощенную модель намагничения: считают
UAV.
(1)
горные породы и руды «идеальным ферромагнетиком». Это дает возможность связать намагниченность с магнитными свойствами. Традиционно прямая задача магнитной съемки заключается в нахождении внешних и внутренних потенциалов для различных объектов. Эти потенциалы удовлетворяют следующему уравнению:
Граничные условия для постановки краевой задачи вытекают из непрерывности потенциала на границе, что вызвано непрерывностью тангенциальной компоненты напряженности, а также из непрерывности нормальных компонент вектора индукции поля.
Таким образом, прямая задача сводится к решению задачи Дирихле для уравнения Пуассона при условии, что все исследуемые объекты находятся внутри расчетной области, и границы этой области достаточно удалены от объектов.
Алгоритм решения прямой задачи сводится к дискретизации уравнения (2) конечно-разностным методом на неравномерной сетке с использованием однородной консервативной схемы [2], [3] и последующей его симметризации [4]. На выходе получается система уравнений вида
где А - действительная, симметричная, сильно разреженная, положительно определенная матрица со слабым преобладанием диагонали. Правая часть уравнения (2) представляет собой дивергенцию остаточной или полной намагниченности.
Алгоритм тестировался на нескольких моделях, имеющих точное решение. Приведем пример расчетов следующей модели: шар радиусом Я = 1 000 м находится в центре некоторой области. Необходимо рассчитать компоненты магнитного поля Вх, Ву, В2 на высоте 50 метров над шаром. Правая часть состоит из компонент: 1х = 1у = 0; Ь = 1 СИ. Для данной модели известны аналитические решения, с которыми и было произведено сравнение модельных. Как видно из рис. 1, рассчитанные и аналитические значения сходятся. Для наглядности выбрана плоскость X = У, в которой Нх и Ну, а также Вх и Ву соответственно, совпадают по абсолютному значению.
При переходе к решению обратной задачи возникает необходимость построения функционала невязки между сигналом наблюденных данных и тестируемой модели с последующей его минимизацией для получения решения.
При решении обратной задачи есть две основные проблемы: высокая расчетная стоимость инверсии и неоднозначность решения. Первая возникает в по двум основным причинам: высокие размерности модельного пространства (томографический подход) и непосредственно высокие расчетные требования прямой 3-0 задачи. Стохастические методы Монте-Карло, позволяющие опре-
сН\gradи=-сНу/ .
(2)
Ах =
(3)
делить апостериорные плотности вероятностей модельных параметров, требуют крайне высокого числа решений прямой задачи даже с малым количеством подбираемых модельных параметров, что делает их неприемлемыми в данной задаче. Вторая проблема связана с минимумом или отсутствием априорной информации и, зачастую, малым количеством данных измерений. Она нивелируется за счет большого количества измерений с БПЛА с разных высот, а также применением техник регуляризации при минимизационном подходе к решению, а также анализе Гессиана функционала невязки с целью оценки качества определения каждого индивидуального модельного параметра.
800
Рис. 1. Тестирование алгоритма решения прямой задачи магнитной съемки
на примере модели сферического тела
Большинство методик решения задачи нацелены на поиски минимума функционала невязки различными методами: прямой подбор, градиентный спуск, квазиньютоновская минимизация и т. п. Для успешного выполнения данных методик необходимо иметь возможность быстро и качественно рассчитывать градиент функционала невязки в пространстве модельных параметров. В качестве одной из самых эффективных техник быстрого расчета производных функционала невязки выступает метод сопряжения [5, 6]. На рис. 2 представлена блок-схема реализованного программно-алгоритмического обеспечения для задачи беспилотной геомагнитной томографии. Синим выделено основное ядро алгоритма: решатель прямой задачи, преобразователь данных измерений во внутренний формат, расчет функционала невязки с учетом погрешностей индивидуальных измерений, применение методов минимизации и регуляризации, итерационная схема цикла решения. Оранжевым - параметры, идущие на вход
и не меняющиеся в процессе решения задачи - априорная информация, трехмерные карты рельефа, непосредственно измерения и их индивидуальные погрешности.
Априорная ■ Стартовая
геология ■ модель
ЗР-рельеф ■ ■ Прямая задача
Тестовый сигнал
Градиент невязки
Измерения (БПЛА)
Разновысотные данные
^_данн
13КИ |Н|
Функционал невязки
Минимизация и регуляризация
Итерационная модель
Погрешности измерений
Невязка итерации
Требуемая невязка
Решение
Рис. 2. Блок-схема программно-алгоритмического обеспечения инверсии
данных разновысотной магнитной съемки с учетом рельефа:
синим выделены внутренние элементы программы; оранжевым - исходные данные и параметризация задачи
В качестве тестовой модели небольшого разрешения нами было выбрано разбиение нижнего полупространства на 27 блоков-параллелепипедов (3 х 3 х 3) с «шахматным» чередованием магнитных восприимчивостей 4 *1 0 _ 5 и 6 *1 0 _ 5, дневная поверхность - синусоидальная функция с полным периодом 1/10 от общей длины профиля по координате X. Рассчитаны тестовые сигналы на трех высотах - 4, 6 и 8 м; к ним добавлен гауссов шум, основанный на модели погрешностей (погрешность увеличивается с высотой и влияет на стандартное отклонение зашумления).
При томографическом подходе геометрия среды является зафиксированной, и задача состоит в восстановлении распределения изучаемых модельных параметров на заранее заданной сетке. В роли стартовой модели выступила однородная среда с постоянным значением магнитной восприимчивости . Для расчета градиента функционала невязки применялся конечно-разностный подход, и тестировался метод сопряженного оператора. В качестве методов минимизации использовались Монте-Карло, Нелдер-Мид, градиентный спуск и квазиньютоновский подход Бройдена-Флетчера-Голдфарба-Шанно (BFGS) в условиях ограниченной памяти (L-BFGS) [7], который и показал наилучший результат минимизации за отведенное время на тестирование. Результаты решения задачи инверсии с использованием L-BFGS представлены на рис. 3. Вы-
несены параметры тестовой модели, итерационная модель инверсии после ограниченного количества шагов (в данном случае 50, время расчета ~2 ч). Выполнена оценка погрешностей: для текущей итерационной модели общая среднеквадратичная невязка сигнала от модели составляет 2,31 %, а отклонение индивидуальных параметров не превышает 5,5 %, что является достаточным для признания модели решением обратной задачи в текущих условиях. Кроме базовых тестов для модели класса «шахматы» проведены эксперименты по определению референтной среды (1 модельный параметр - полупространство) и оценка влияния рельефа на получаемый сигнал измерений.
6,5
6,15
6,27 6,29
5
■
о
1
X
5,97 5,98 5,99 6,04 6
6 { 5*5 5,84 5,82 5,.89 •
5,75
6,02 5,89 5,9
5,5 5 4,5 4 3,5
3,98 3,96 3,96 • • •
4 2 4 21
4,1 4,2 4,21 4,11
3,92
I
4,1
3,87
• Исходная модель
• Инверсия
4,12
3,78 3,76 • • •
3
10 13 16 19 Номер параметра модели, п
22
25
Рис. 3. Тестируемая модель, итерационная модель. Сравнение
модельных параметров
1
4
7
Полученные результаты подтверждают, что предложенная схема инверсии работает для синтетической модели небольшого разрешения и открывает возможность перейти к тестам данных реальных измерений и большего числа модельных параметров, анализу матрицы чувствительности данных к модельным параметрам; где в последствии добавятся реальные погрешности измерений, неоднородность сетки наблюдения, априорная информация о геологии.
Работа выполнена в рамках проекта «Разработка программно-алгоритмического обеспечения для технологии беспилотной геомагнитной томографии в условиях криолитозоны» по программе фундаментальных исследований президиума РАН «Арктика - научные основы новых технологий освоения, сохранения и развития».
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Блох Ю. И. Интерпретация гравитационных и магнитных аномалий. - 2009. - 231 с.
2. Самарский А. А. Теория разностных схем. - М. : Наука, 1979. - 655 с.
3. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. - М. : Наука, 1978. - 592 с.
4. Кузнецов Ю. И., Агапитова Н. С. Математические основы моделирования на ЭВМ. -Ю.-Сахалинск : Из-во ЮСИЭПИ, 2003. - 214 с.
5. Maksimov M., Velimsky J. Fast calculations of the gradient and the Hessian in the timedomain global electromagnetic induction inverse problem // Geophys J Int (2017) - 210 (1) -P. 270-283.
6. Fletcher R. Practical methods of optimization (2nd ed.) - New York : John Wiley & Sons,
1987.
7. Fichtner A. Full Seismic Waveform Modelling and Inversion. - AGEM, 2011.
REFERENCES
1. Bloh Yu. I. Interpretaciya gravitacionnyh i magnitnyh anomalij. - 2009. - 231 s.
2. Samarskij A. A. Teoriya raznostnyh skhem. - M. : Nauka, 1979. - 655 s.
3. Samarskij A. A., Nikolaev E. S. Metody resheniya setochnyh uravnenij. - M. : Nauka, 1978. - 592 s.
4. Kuznecov Yu. I., Agapitova N. S. Matematicheskie osnovy modelirovaniya na EVM. -Yu.-Sahalinsk : Iz-vo YuSIEPI, 2003. - 214 c.
5. Maksimov M., Velimsky J. Fast calculations of the gradient and the Hessian in the timedomain global electromagnetic induction inverse problem // Geophys J Int (2017) - 210 (1) -P. 270-283.
6. Fletcher R. Practical methods of optimization (2nd ed.) - New York: John Wiley & Sons,
1987.
7. Fichtner A. Full Seismic Waveform Modelling and Inversion. - AGEM, 2011.
© М. А. Максимов, И. В. Суродина, В. Н. Глинских, 2018
