Бесплатный автореферат и диссертация по геологии на тему

Алгоритмы расчета аномального поля 2D и 3D сильно намагниченных тел и их реализация в среде Unix

ВАК РФ 04.00.12, Геофизические методы поисков и разведки месторождений полезных ископаемых

Текст научной работыДиссертация по геологии, кандидата технических наук, Костров, Николай Павлович, Екатеринбург

/

Российская академия наук Уральское отделение Институт геофизики

Костров Николай Павлович

Алгоритмы расчета аномального поля 2D и 3D сильно намагниченных тел и их реализация в среде Unix

Специальность 04.00.12 - Геофизические методы поисков и разведка месторождений полезных ископаемых

Диссертация на соискание ученой степени кандидата технических наук

Научные руководители:

доктор геолого-минералогических наук профессор В.В. Кормильцев, кандидат физико-математических наук В. А. Шапиро.

Екатеринбург 1999

Оглавление.............................................................................................2

Введение....................................................................................................3

Глава 1. Алгоритм определения магнитного поля неоднородных тел, основанный

на объемном векторном интегральном уравнении......................................6

1.1. Объемные и поверхностные интегральные уравнения, особенности подходов при численном решении..............................................»....„„б

1.2. Вывод интегрального уравнения для напряженности магнитного поля магнетика в двумерном случае. Трехмерный аналог...............................13

1.3. Дискретизация интегрального уравнения, связь с теоремой Пуассона........15

1.4. Компоненты тензора Грина в двумерном случае, вычисление в комплексной области......................................................................18

Глава 2. Исследование алгоритма....................................................................22

2.1. Аналитические модели для тестирования алгоритма..............................22

2.2. Анализ кривых насыщения в двумерном случае.............................. ,,...26

2.3. Анализ кривых насыщения в трехмерном случае.................................29

2.4. Анализ симметрии внутреннего поля.................................................31

2.5. Исследование устойчивости внутреннего и внешнего поля при

большом числе разбиений....................................................................32

2.6. Разрушение внутреннего поля.........................................................34

Глава 3. Программные комплексы для 2D и 3D случаев в ОС Unix.........................38

3.1. Операционная система Unix.............................................................39

3.2. Пакет программ МАГЛАБ-П для двумерного моделирования,,,...............40

3.3. Некоторые алгоритмы вычислительной геометрии, использованные

при создании комплекса Маглаб-П........................................................45

3.4. Моделирование аномального поля в горной выработке.........................49

3.5. Примеры подбора наблюденного поля.............................................53

3.6. Пакет программ МАГЛАБ для трехмерного моделирования...................58

3.7. Моделирование эффекта подмагничивания Манчажской региональной магнитной аномалии вариациями земного поля................61

Заключение..............................................................................................72

Литеретура..............................................................................................72

Введение.

Актуальность исследования

Уровень исследований современной магнитометрии непрерывно улучшается, появляются новая высокочувствительная и стабильная аппаратура, новые методы навигации, методики учета естественных и искусственных помех, что позволяет получать значения аномального магнитного поля с недоступной ранее точностью. Неотъемлемой частью магниторазведки в наши дни является широкое применение современных машинных методов решения прямых и обратных задач. Необходимо дальнейшее развитие алгоритмического и программного обеспечения для численного решения прямых задач магнитометрии и интерпретации методом подбора в случае неоднородной магнитной восприимчивости рудных и некоторых региональных объектов. Необходимо внедрение в практику математического моделирования свободно распространяемых операционных систем, альтернативных MS DOS и Windows, имеющих некоммерческий характер и создающих обширные возможности для научно-технических расчетов.

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

Задачи исследования.

1. Разработка нового алгоритма для вычисления магнитного поля от совокупности 2D многоугольных неоднородных по магнитной восприимчивости призм на основе объемного векторного интегрального уравнения для напряженности.

2. Создание программных комплексов в ОС U№X для моделирования напряженности магнитного поля неоднородных 2D и 3D объектов.

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

4. Иллюстрация эффективности разработанного алгоритма и пакетов программ на практических примерах.

Научная новизна.

1. Для интегрального уравнения Фредгольма 2-го рода, описывающего напряженность сильно намагниченного тела или совокупности тел с неоднородными магнитными свойствами в двумерном случае, выполнена дискретизация и предложен алгоритм вычисления компонентов тензора Грина для треугольного сечения, сочетающий интегрирование в плоскости действительных переменных с формулой A.B. Цирульского для внешнего логарифмического потенциала многоугольного контура.

2. Предложенный алгоритм реализован в системе научных расчетов SCILAB в ОС UNIX, основанной на принципах модульности и открытости, что позволяет рационально организовать взаимодействие с программным комплексом, и осуществлять видоизменение и развитие программ.

Основные защищаемые положения сформулированы в заключении диссертации.

Практическая значимость.

Разработанные пакеты программ Маглаб и Маглаб-П позволяют получить численные решения для напряженности магнитного поля совокупности неоднородно и сильно намагниченных 2Е) и ЗБ объектов. Пакеты могут быть применены для изучения многих проблем магнитометрических исследований. В качестве примера приведем следующие:

1. Исследование практической эквивалентности аномального поля над объектами сложного внутреннего строения.

2. Исследование вопросов подмагничивания на магнитных аномалиях для оценки <3 фактора. В этом случае ЗБ моделирование является неотъемлемым технологическим звеном при планировании и оценке результатов эксперимента в реальном масштабе времени.

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

Реализация результатов работы.

В настоящее время пакеты программ применяются в исследовательских целях в Институте геофизики УрО РАН в проекте РФФИ N98-05-64816 "Комплексное геофизическое и геологическое изучение структуры зон сочленения

палеоконтинентальных и палеоостроводужных террейнов на примере Южного Урала" для переинтерпретации результатов магнитометрии в Сакмарской структурно-тектонической зоне и на Главном Уральском Глубинном Разломе, а также в Уральской государственной горно-геологической академии при подготовке квалификационных работ инженеров и магистров, в РФЯЦ - ВНИИТФ(г. Снежинск) при изучении проблемы использования магнитометрии для поисков и локализации подземных техногенных явлений.

Апробация работы.

Результаты исследований и основные положения диссертации докладывались на 8-th Scientific Assembly of IAGA with ICMA and STP Symposia. 4-15 Aug. 1997,Uppsala, Sweden, на 67-th Annual Meeting of Society of Exploration Geophysicists, New Orleans, Louisiana, September 13-18, 1998, на международной конференции "Проблемы геодинамики, сейсмичности и минерагении подвижных поясов и платформенных областей литосферы", Екатеринбург, сентябрь-октябрь, 1998, на Международном семинаре им. Д.Г.Успенского "Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей" (26 сессия, Екатеринбург,25 -30 января 1999 г.). По теме диссертации опубликовано 9 статей и 8 тезисов докладов.

Личный вклад автора.

Работа выполнена в лабораториях магнитометрии и электрометрии Института геофизики УрО РАН в 1994-1999 г.г. Алгоритмические основы программ и приложения разработаны совместно с научными руководителями, программный продукт - лично автором диссертации. Программные продукты(пакеты программ Маглаб и Маглаб-П) предоставляются бесплатно всем пользователям.

Объем работы.

Диссертация изложена на 78 страницах с 29 рисунками и состоит из введения, трех глав, заключения и списка литературы из 72 названий.

Автор признателен научным руководителям д.г.-м.н. профессору В.В. Кормильцеву и к.ф.-м.н. В. А. Шапиро, своему соавтору по основным публикациям к.т.н. А.Н. Ратушняку, а также к.ф-м.н. Н.В. Федоровой, д.ф.-м.н. профессору П.С. Мартышко, к.г.-м.н. А.Н. Бахвалову, к.г.-м.н. A.B. Чурсину, к.ф.-м.н. Ю.К. Доломанскому за

консультации, предоставление материалов, обсуждение и рецензирование результатов и м.н.с. Н.Н. Винничук за помощь в практических расчетах.

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

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

1.1. Объемные и поверхностные интегральные уравнения, особенности подходов при численном решении.

Существует два основных подхода к численному моделированию. На основе дифференциальных уравнений разработан метод сеток, для реализации которого в модель должна быть включена значительная часть пространства, окружающего неоднородности. Преимущество интегральных уравнений в том, что искомые поля вычисляют интегрированием только по аномальным областям, а не по всему пространству. Методы численного решения интегральных уравнений Фредгольма 2-го рода в геофизике, главным образом, развивались при решении задач вызванной поляризации и расчете электромагнитного поля 3-D тел. Интегральное уравнение получали путем записи системы уравнений Максвелла [1] для падающего и рассеянного электромагнитного поля и решения их методом функции Грина [2,3]. Компоненты Gу (г,г') тензора Грина выражают /-ый компонент электрического поля в точке Р(г),

который возбуждается элементарным электрическим диполем, расположенным в точке б(г') внутри объема V и ориентированным вдоль оси х ■ прямоугольной системы

координат [4]. Согласно работам [4,5,6,7,8], интегральное уравнение для электрического поля Е(г) записывали как:

Е(г) = Е,(г) + (а-2-о-1)|С(г,г')Е(г')^Г, (1.1.1)

V

где ЕДг)- падающее электрическое поле, 6'(г,г')- функция Грина, сГрО^-

электропроводность среды и неоднородности соответственно. Следует отметить, что в этом уравнении пренебрегают неоднородностью свойств аномального объема V как и в работах [6,8,11,12,13,14,16,17]. В работах же [4,7,10,15,18-24] такую неоднородность учитывают, что соответствует внесению множителя - под знак интеграла. Далее

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

Е(г) = Е,(г) + (<т2-о-р^Ея |(7(г,г') йУ\ (1.1.2)

/»=1 ¥п

- г -ИсхК _

Путем дальнейших преобразований уравнение (1.1.2) сводилось к системе линейных алгебраических уравнений (СЛАУ) для вектора напряженности электрического поля в элементарных ячейках:

ТУ

I

п=1

сг0 -<7, I 1Г _х

1 тп °тп

°1

Ей =-Еда. (1-1.3)

где 8тп - символ Кронекера, Ттп - интеграл от функции Грина. После того, как решение СЛАУ (1.1.3) получено, можно вычислить Е(г) вне тела по формуле (1.1.1). Таким образом, в данном подходе основную проблему представляет вычисление интеграла от функции Грина. В работе [7] интегралы от недиагональных элементов тензора Грина вычисляли методом центральной аппроксимации, суть которого в следующей замене:

0(г,г')й?Г -» С(г,г') \<ЯГ . (1.1.4)

п "п

При достаточно больших Ыэтот способ дает адекватные результаты [7]. В работе [8] найдено, что для ячейки в форме прямоугольного параллелепипеда такой метод применим для расстояний между центром текущей ячейки и центрами остальных ячеек Кпт>5-Мт(Ах,Ау,Аг), (1,1.5)

где Лх,Ау,Лл - размеры ячейки по осям координат. Для меньших расстояний предлагалось проводить численное интегрирование. В работе [7] элементарные ячейки Уп заменялись сферами равного объема, что давало возможность аналитически

вычислить значение тензора Грина для диагональных элементов, но существенное отличие формы ячейки от сферы приводило к плохим численным результатам. Поэтому предлагалось численно интегрировать функцию Грина по дополнению сферы до ячейки и добавлять результат к аналитическому результату для сферы. В работе [8] вычисления интеграла от функции Грина предлагалось выполнять численно. В работе [6] брали кубические ячейки и функцию Грина разделяли на части для исключения проблем со сходимостью из-за двойных производных (2) и учета влияния границы земля-воздух. Выделяли часть, возникающую из-за токов С}^ и часть, соответствующую

положительным и отрицательным зарядам Ср. Распределение зарядов имеет место на

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

тогда как вторичная часть (индекс Б) учитывает границу земля-воздух:

СА=0А+0А> °<р = °!р+(}р- ■ (1Л-6>

р

Кубическая ячейка заменялась сферической равного объема и интеграл от функции Сд

вычислялся аналитически отдельно для диагональных и недиагональных компонентов. Для вычисления О у интегрирование по объему элементарной ячейки заменяли на

интегрирование по поверхности объема. Для аналитического вычисления диагональных Р

элементов квадратную грань аппроксимировали круглым диском той же площади.

р

Для вычисления недиагональных элементов О^ применяли разложение

Цд+1)

1 .

-в ряд -+—---1——, что соответствует низким частотам и

4 я 1С 4я1С \Ъс

о

коротким расстояниям. Части функции Грина, учитывающие границу земля-воздух Од и

о

Сф, вычисляли методом численного интегрирования с применением функций Бесселя 0-го и 1-го порядка [9]. В работе [10] авторы модифицировали схему расчета, заменив

концентрацию заряда на границах между ячейками на однородное объемное распределение заряда от центра одной ячейки до центра следующей.

Исходной при решении магнитостатической задачи с учетом размагничивания является следующая краевая задача[11]. Пусть на границе дО односвязной области Г) евклидова пространства Яп выполняются следующие условия

Д^=0, к2,

У1=У2, (1.1.7)

д¥? дУ, дЖ

г дп 1 дп 1 А дп где У~2 и Уу внутренний и внешний потенциал области О, /л^ и /л^ - магнитная

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

В работах [И, 12,13] предложено решение методом линейного сопряжения прямой задачи магниторазведки с учетом размагничения. Суть метода в следующем. Перейдем к комплексным переменным ¿=х + гу, 2-х-1у и определим следующие фукнции

Мг) = -—

К п

Гдук _

дх ду

-I-

дх ду

к = 1,2, и^г) = -

Тогда задача (7) сводится к задаче линейного сопряжения [13] относительно функций

(1.1.8)

яс^до, (1.1.9)

¿и7 + щ

где /'(О - производная правой части уравнения кривой дВ, Х-———Далее

решение искали в форме = —— [ ^^ ¿/^. В результате некоторых

преобразований было получено интегральное уравнение для искомой плотности интеграла типа Коши Ч'(г) :

ш П-Ч + АУ) Ат 1 г ЧЧО ^ "3/~"3

у = ^—-¿ +(р> АЧ> = — — (р = —-——. (1.1.10)

я ^аЬ^о л

Это уравнение решали методом последовательных приближений. В [11] также указан способ решения с обходом особой точки ^.

В работе [14] предложено решение прямой трехмерной задачи магниторазведки с учетом размагничения путем сведения задачи (1.1.7) к интегральному уравнению относительно искомой плотности двойного слоя а:

К.п 2яЛ у <"1 т) I .! /1 1 ил

=--, К = г— г*, (1.1.11)

2 п /¿2 + М\

<г(т)~\<т(г')

¿71 8

ЯЪ 5

где п- внешняя нормаль к поверхности ¿'.Решение находилось методом последовательных приближений и по найденной плотности а вычислялось аномальное магнитное поле:

Н = (1.1.12)

5 К

В работе [11] на основе уравнения (1.1.11) получено для класса тел, звездных относительно некоторой внутренней точки, оптимальное решение данной задачи: вычисление поверхностного интеграла было