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

Кокшаров, Дмитрий Евгеньевич

Бесплатный автореферат и диссертация по наукам о земле на тему
Алгоритмы и новые компьютерные технологии решения структурных обратных задач гравиметрии и магнитометрии
ВАК РФ 25.00.10, Геофизика, геофизические методы поисков полезных ископаемых

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

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

Кокшаров Дмитрий Евгеньевич

АЛГОРИТМЫ И НОВЫЕ КОМПЬЮТЕРНЫЕ ТЕХНОЛОГИИ РЕШЕНИЯ СТРУКТУРНЫХ ОБРАТНЫХ ЗАДАЧ ГРАВИМЕТРИИ И МАГНИТОМЕТРИИ

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

25.00.10-Геофизика, геофизические методы поисков полезных ископаемых

Екатеринбург - 2005

Работа выполнена в Институте геофизики УрО РАН

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

доктор физико-математических наук, профессор Мартышко Петр Сергеевич

Официальные оппоненты:

доктор физико-математических наук Романюк Татьяна Валентиновна

доктор геолого-минералогических наук Иголкина Галина Валентиновна

Ведущая организация:

Уральский государственный горный университет, г. Екатеринбург.

Защита состоится 23 декабря 2005 г. в 14.00 на заседании диссертационного совета Д 004.009.01 при Институте геофизики УрО РАН, по адресу: 620016, г. Екатеринбург, ул. Амундсена, д. 100.

С диссертацией можно ознакомиться в библиотеке Института геофизики УрО РАН.

Автореферат разослан «¿£>> ноября 2005 г.

Ученый секретарь диссертационного со доктор физико-математических наук

1247311

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы

Важнейшими способами исследования структуры недр Земли в геофизике являются методы грави- и магниторазведки. С их помощью ведется изучение плотности геологических пород и их магнитной проницаемости. Измерения аномалий гравитационного и магнитного поля Земли используются при поиске и разведке месторождений полезных ископаемых, для изучения приповерхностных структур Земли, а также для решения экологических задач. В теории потенциальных геофизических полей различают два важнейших класса задач - прямые и обратные. Прямая задача - нахождение поля, порождаемого геологической структурой с известными свойствами (параметрами): распределением плотности или намагниченности. Такие задачи (моделирования полей) достаточно хорошо исследованы как для гравитационных, так и для магнитных полей в двух- и трехмерном случае (можно отметить работы В.Н. Страхова, Ю.И. Блоха, П.С. Мартышко, С.М. Оганесяна, В.И. Старостенко, а также М. Talwani, S.P. Sharma, Р. Weidelt и других). Обратные задачи - определение физических параметров пород по измеренным полям. При разработке теории решения обратных задач гравиметрии и магнитометрии лучше всего исследован двумерный случай. Результаты В.Н. Страхова, A.B. Цирульского, В.К. Иванова, П.С. Мартышко, Н.В. Федоровой, Ф.И. Никоновой позволяют эффективно находить решение задачи - эквивалентное семейство решений (существенный прогресс в этом направлении достигнут с применением результатов ТФКП). В трехмерном случае обратные задачи решаются менее успешно. Причин этому несколько. Во-первых, такие задачи сводятся к решению обратных задач математической физики, а последние, в свою очередь, являются некорректно поставленными, что существенно затрудняет интерпретацию полей. Общая теория и методы решения некорректно поставленных задач были разработаны А.Н. Тихоновым,

B.K. Ивановым, M.M. Лаврентьевым, а применительно к задачам геофизики методы регуляризации разрабатывали В.Н. Страхов, В.Б. Гласко. Большой вклад в разработку теории и практическое применение методов решения трехмерных обратных задач внесли А.И. Кобрунов, В.М. Новоселицкий, В.Н. Страхов, В.И. Старостенко, П.С. Мартышко, И.Л. Пруткин, A.C. Долгаль. Во-вторых (это следует из результатов по теории эквивалентности В.Н. Страхова, A.B. Цирульского, С.М. Оганесяна, В.И. Старостенко и других), решить однозначно обратную задачу при интерпретации геофизических полей естественной природы без дополнительной информации в принципе невозможно. Кроме того, задачи такого класса требуют, как правило, существенных вычислительных ресурсов, которые до последнего времени были недоступны. Следовательно, решение задач практической интерпретации данных потенциальных геофизических полей остаются актуальными, особенно с привлечением современных вычислительных средств. Связь работы с научными программами, планами, темами. Работа выполнялась в соответствии с плановой тематикой НИР Института геофизики УрО РАН по теме 01.200.2 09053.

Цель и задачи исследования.

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

Предмет и объект исследования.

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

Научная новизна полученных результатов.

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

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

В рамках одного пакета осуществляется полный цикл обработки данных:

• загрузка исходных файлов;

• исключение боковых (находящихся вне рассматриваемой области) источников поля;

• выделение поля от аномалий, залегающих в изучаемом слое;

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

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

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

Практическая значимость полученных результатов.

Компьютерная технология решения линейной трехмерной обратной задачи для плоских слоев и слоев, границы которых являются поверхностями, имеющими горизонтальную плоскую асимптоту, была успешно применена для интерпретации практических данных по 3 регионам, результаты были переданы заказчикам и включены в отчеты ЗАО «Пермь-Лукойл» (г. Пермь) и ЗАО «Гравиразведка» (г. Москва).

Сведения об апробации результатов диссертации.

Результаты научной работы доложены на российских и международных

конференциях:

1. 66th EAGE Conference & Exhibition. Extended abstracts. Paris, 2004.

2. VI междисциплинарного международного симпозиума «Строение и эволюция геосфер». Хабаровск, 2003.

3. 36 Региональная молодежная конференция «Проблемы теоретической и прикладной математики», ИММ УрО РАН, 2005.

4. Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей. Международный семинар им Д.Г. Успенского, 2004.

5. 6 Уральская молодежная научная школа по геофизике. Горный институт УрО РАН. Пермь, 2005.

6. Глубинное строение, геодинамика, мониторинг, тепловое поле Земли, Интерпретация геофизических полей. Третьи научные чтения

Ю.П. Булашевича. 2005.

7. VII международный симпозиум «Закономерности строения и эволюции геосфер», Владивосток, 2005.

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

1. Для трехмерных задач гравиметрии разработаны новые алгоритмы решения в линейной и нелинейной постановке:

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

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

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

2. Интегродифференциальное уравнение теоретической обратной задачи о магнитных контактах с учетом размагничивания сведено к системе нелинейных алгебраических уравнений.

Структура и объем диссертации: диссертация состоит введения, 5 глав, заключения и списка литературы. Общий объем работы составляет 85 страниц, включая 37 рисунков, библиографический список содержит 103 наименования.

Работа частично (компьютерная технология в методе поиска структурных границ) поддержана грантом УрО РАН для молодых ученых за 2005 г. Автор выражает глубокую благодарность научному руководителю - д.ф.-м.н., профессору П.С. Мартышко, а также д.ф.-м.н. И.Л. Пруткину, чл.-корр. РАН, профессору В.В. Васину, к.ф.-м.н. E.H. Акимовой, к.ф.-м.н. Н.В. Федоровой, к.ф.-м.н. Ф.И. Никоновой за консультации и интересные обсуждения в процессе работы.

СОДЕРЖАНИЕ ДИССЕРТАЦИОННОЙ РАБОТЫ

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

В главе 1 (Методы решения обратных задач гравиметрии) кратко излагаются современные подходы к математической интерпретации гравитационных аномалий.

1. Интерпретация гравитационной аномалии методом подбора (на основе решения прямой задачи).

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

3. Подбор измеренных данных полем от сингулярных источников.

4. Методы инверсии оператора решения прямой задачи.

Глава 2. Выделение аномалий поля в указанной области, подготовка данных

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

Подготовка данных

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

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

где хьуь - к-я точка, в которой определено значение Требуется найти функцию ^х,у), такую, чтобы норма разности ||^хк,ук)-Гк|| была мала. При этом на интерполирующую функцию фс,у) могут налагаться некоторые дополнительные ограничения (как правило, требования гладкости или ограничения на скорость роста). В работе изложены основы применявшихся в процессе подготовки практических данных методов: Шепарда, разложения по базисным функциям (мультиквадрик) триангуляции, кригинга, метода наименьшей кривизны. Как показали эксперименты с реальными данными измерений геофизических полей, технология предварительной подготовки данных играет важную роль в процессе интерпретации, в частности, неаккуратный выбор метод интерполяции может приводить к созданию мнимых экстремумов функции, которые будут впоследствии перенесены в этап

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

Выделение аномалий поля в указанной области

Как показывают исследования геологических структур, существенный вклад в гравитационное поле исследуемого слоя могут оказывать источники, залегающие вне него. В частности, это источники, находящиеся вне рассматриваемой области в проекции на земную поверхность (боковые источники), а также источники, залегающие на глубинах меньше и больше рассматриваемой. Для решения задачи об определении плотности в слое необходимо выделить поле от этого слоя, производя следующие операции (И.Л. Пруткин, П.С. Мартышко, 1998).

1. Исключение боковых источников поля. Производится решение двумерной задачи Дирихле (краевыми условиями для задачи Дирихле являются значения поля на краях исследуемой области) Ди(х,у) = 0

и(х,у)|ао=АЕю"(х,у)|лз (21)

Затем найденное решение задачи Дирихле вычитается из исходных значений поля.

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

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

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

/(Л

4Ы

г /-.чЛ

сЬсф'

что позволяет говорить о том, что гармоническая функция отражает плавные изменения регионального фона внутри области.

При численной реализации метода в программе поиск гармонической функции и осуществлялся решением методом сеток дискрет из ированной системы (2.1), оператор Лапласа заменялся сеточным оператором

где Ь=Дх=Ду, либо

Аи(х,,у )= '+ """ + ' -О,

J Ах Ду

при различии шагов по х, у олученная система дополнялась уравнениями краевых условий що=А£т"(х,у),..., Дё™" (х, у) и решалась относительно иу.

2. Выделение источников гравитационного поля, локализованных в горизонтальном слое между указанными глубинами. Как правило, из априорных данных (например, результатов сейсморазведки) известно, что исследуемые аномалии геологических структур залегают на определенной глубине под дневной поверхностью. Модель латерально неоднородной среды была предложена В.М. Новоселицким. С целью выделения из общего поля составляющей от источников в указанном слое, применяется следующий метод. С помощью формулы Пуассона

и(х',у',Н) = ~ Г Г-------^«(^,0)^,(2.2)

дающей решение краевой задачи Дирихле

Ди = 0,г>0 и = и(х,у,0),2 = 0

для полупространства, производят пересчет величины поля на высоты Н( и Н2, затем поле с помощью этой формулы продолжается на глубину -Н; (в формуле 2.2 производится замена Н на -2Щ с регуляризацией (А+аЕ). Затем поле пересчитывается по формуле (2.2) на земную поверхность. При пересчете вниз по (2.2) существенно уменьшается влияние источников, залегающих в слое ге[0,-Н;], поскольку (2.2) предполагает отсутствие источников поля в области применения (при наличии аномалий плотности в них нарушается условие Аи = 0).

Затем находится разность между двумя вычисленными значениями поля. Это разностное поле принимается за поле от источников внутри указанного слоя и по нему производится решение обратной задачи в какой-либо постановке. Можно привести грубую оценку ослабления поля от источников, которое имеет место при пересчете поля вверх по формуле (2.2). Рассмотрим изменение максимума аномального поля от одиночного шара радиуса г с избыточной плотностью До, центр которого залегает на глубине Н, при перечете поля на

47гг3Д(Т 4яг3Да высоту ДН (Г - гравитационная постоянная): к = 5-^-.

ЗН2 3(Н + ДН)2

откуда при ДН = 2Н (шар, залегающий в середине исключаемого слоя) ослабление составит к=9, что довольно существенно.

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

например, для сфероидов гравитационный эффект представляется формулой (Голов И.Н., Сизиков B.C., 2005):

где у - гравитационная постоянная, е=(1-(с/а)2)ш - эксцентриситет, Д<7 -избыточная плотность, (х-о,Уо,2о) - координаты центра сфероида, а=Ь -горизонтальные, с - вертикальная полуось сфероида, причем а < с.

Если теперь рассмотреть модель, в которой сфероид со сферичностью В = —

(е=1 - шар) расположен в слое толщины 4с, так, что его центр находится на глубине 2с (посередине слоя), и найти зависимость ослабления поля Г при

пересчете на высоту 4с (толщину слоя) от сжатия эллипсоида Б = —,

то получится следующий график (рис. 2.1).

&g(x>y)=^щАа-^jip-arctgp )z(

ае

■о г

г 4 б I г и

Рис. 2.1.

Из графика (рис. 2.1) хорошо видно, что уже при е = 3 ослабление принимает величину порядка 5 (что еще допустимо), а при £ = 5 приблизительно равно 3, что практически неотличимо от «нормального» поля «глубоких» источников. Это означает, что вытянутые в горизонтальном направлении объекты (эмпирически подобранное соотношение высоты к горизонтальному размеру более 1:3) слабо отсеиваются методом пересчета. Для исключения таких объектов необходимо использовать дополнительные методы -предварительный подбор сингулярными источниками (например, стержнями), с последующим вычитанием вычисленного поля подобранной модели из измеренного поля.

Глава 3. Обратная задача гравиметрии

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

Структурные границы. Метод локальных поправок

Модель этого варианта геологической структуры - структурной границы (или

«контактной поверхности») описывается следующими формулами. Под

дневной поверхностью между глубинами h- и h+ залегает поверхность S,

разделяющая два слоя плотности а, и сг2 соответственно. Поверхность S

задается уравнением z=z(x,y) (требуется однозначность, и, по необходимости,

достаточная гладкость функции z(x,y)). При этом выполняется предел

,'jm \2(х> у)~н\- О,

М-»»'

М-»»

то есть у поверхности S существует горизонтальная асимптота z=H (очевидно, h_<H<h+).

Поле структурной границы (контактной поверхности) описывается следующей формулой:

Дг(*',/,0) = Д<7/| \

1 1

= ¿хЛу,

—ос —ОС

^х-х')2+(у-у')2 + г(х,у)2 ^(х-х')2+(у-у')2+И\

где Д£(х,у,0) - вычисленное разностное поле на дневной поверхности, Да=СТ2-сг 1 — разность плотностей между слоями, Г - гравитационная постоянная, И -глубина залегания асимптоты структурной границы, г(х,у) - глубина залегания собственно структурной границы. Эта формула описывает разность между полем от асимптоты и от границы раздела. Она непосредственно следует из

внести производную под знак интеграла и вычислить внутренний интеграл.

Существенно упрощает ситуацию тот факт, что гравитационное поле обладает свойством аддитивности, это позволяет, в частности, порознь считать поля от различных аномалий, рассматривать множественные структурные границы. Одним из недавно разработанных методов решения задачи о структурных границах (Пруткин И.Л., 1986; Мартышко П.С., Пруткин И.Л., 2003) является метод локальных поправок. Основой этого метода является предположение о том, что на значение поля в данной точке дневной поверхности наибольшее влияние оказывает значение глубины залегания ближайшей точки контактной поверхности.

Формула вычисления поля структурной границы как интегральное уравнение относительно г(х,у):

¿Ьсс1ус1г '

№х',у',0) = А<т/ $ !

Полагаем, что значения поля известны на сетке {хьэд} с равномерным шагом по х и у.

Эту формулу можно рассматривать как нелинейное интегральное уравнение 1-го рода относительно искомой функции г(х,у). Обозначим подынтегральное выражение

К(х', у', х, у,г(х, у)) = , * —, * =--

^(х-х')2+(у-/)2+2(х,у)2 ^(х-х'У+Су-УУ+Л2

Пусть ] - значения искомой функции, полученные на п - ом шаге.

Обозначая

' J

где с - коэффициент кубатурной формулы, г^ = 7(п) (х,., ^), следующее приближение

Тогда (Пруткин И.Л., 1986)

г"** = 1

>! 1+г:г{и„-и;)

Регуляризованная модификация данной формулы имеет вид

г"*1 = 2,1

4 Х + а-грц-и;)

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

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

Слой переменной плотности.

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

А^',/,0) = Ш , 2 tedydz.W

& J(x-x')2+(y-y')2+(z-ï)2

Задача - по измеренному на земной поверхности полю Ag(x',y',0) найти о(х,у). Таким образом, относительно неизвестной функции двух переменных записано линейное интегральное уравнение первого рода. Поскольку аналитическое решение такого уравнения затруднительно, естественно применить дискретизацию. При этом возможны два подхода.

Первый подход предполагает представление слоя в виде набора однородных прямоугольных призм с горизонтальными размерами Дх, Лу и простирающихся по вертикали от Н, до Н2. Призмы однородны, а поле такой призмы единичной плотности Agij(x',y',0) можно точно вычислить по формулам Тальвани. Следовательно, можно записать:

Получается система линейных алгебраических уравнений относительно ст(х,у). Однако, из-за высокой вычислительной сложности формул Тальвани, такой подход потребует для решения больших затрат машинного времени. Второй подход основан на том, что можно свести трехмерный интеграл (3.1) к двумерному следующим образом:

^ Ж во И 2

а?

сг(х,у)

* ' ^(х-х')г+(у-у')2+(2-2'У

¿Ыскуйг •

00 00 Л 2 а

¡д

00 X

-/И

1 1

с1хс!у =

а-{х,у)(Ыу

\л!(х-х')г +(у-у')г +Н* ^(х-х')гНу-у')2+Н2г ,

Значит, справедливо следующее выражение для гравитационного поля от горизонтального слоя переменной плотности:

_1_

д/(*-дОг+(у-у')г+Н1г (3.2)

сг(дг ,у)сЬсс1у

во 00

Ь8(х\у,,Щ = Г\\

^х-х')г+{у-у')г+Нгг )

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

Уравнение (3.2) можно записать в операторном виде Аа = Д£. В этом случае

двумерный массив значений плотности а(х„у,) на сетке 1=1,...,N. ]=1,...,М рассматривается в виде одномерного вектора (например, записывается

построчно) и умножается на матрицу интегрирования размерности МИ х \ДО. Можно показать (Пруткин И.Л., 1988), что оператор уравнения является самосопряженным и положительно определенным, следовательно, для регуляризации возможно применение метода Лаврентьева, то есть вместо решения Асг — д^ искать решение (A + aE)o^ = Ag, где а - параметр регуляризации. Легко видеть, что для такого уравнения дискретизация также приведет к системе линейных алгебраических уравнений. Видно, что именно в этом методе особо важен процесс выделения источников поля, залегающих в указанном слое.

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

На самом деле, можно существенно расширить класс изучаемых функций следующим приемом. Как правило, в случае достаточно однородных пород априори известен закон изменения плотности с глубиной, то есть возможна запись а(х,у,г)=а(х,у)(р(г). Тогда при подстановке в (3.1) имеем:

^ оо оо Н 2

<т(х,у,г)

00 00 л #2

1я, л/(*-*')2 + (у -у')2 +(* -*')' ст{х,у,2)

=сЫусЬ =

со до И 2 л

-ОО-ЖЯ,

д_ дг'

°(х,у)<р(г)

,у!(Х-х')2 + (У - у')2 + (2- 2'У ,

сШуйг =

-//И*.»

1 Я,

(з(г)

¡¿г

То есть формула (3.1) приобретает вид:

¡¿и/у.

оо ао Л} ¿ч

Ф)

В общем случае, интеграл

/(*,**•,/,*•)="} й *«> X

может оказаться не берущимся аналитически, однако, например, в случае ^

функции А^г') полиномиального вида (а с помощью такой функции можно сколь угодно точно приблизить любую непрерывную функцию) такой интеграл ^

можно вычислить явно. После нахождения функции 1(х,у,х\у',2') дальнейшая процедура в точности совпадает с описанной для случая сг=<т(х,у), то есть задача приближенного вычисления сводится к решению системы линейных алгебраических уравнений (в том числе и в случае с регуляризацией).

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

Подобные рассуждения можно провести и для случая, когда граница слоя -произвольная гладкая поверхность. Соответствующая модель предполагает отсутствие аномалий плотности вне слоя трехмерных слоев, границы которых являются поверхностями, имеющими горизонтальную плоскую асимптоту (слоя с «криволинейными» границами) Н^Н^у) и Н2=Н2(х,у), такими, что

Н,>Н2 V(x,y), и выполняется Н,(х,у)—х_>±зд >h, = const, при этом

у-»±<ю

задается распределение плотностей ст=а(х,у) внутри слоя. Тогда (аналогично выводу (3.2))

Ag(x',/,0) = /JJ

Jix-xY+iy-y'y+H^y)2 (3.3) a(x,y)dxcfy

J(x-x')2 +(у-уУ +Н2(х,у)\

Решение обратной задачи в этом случае производится следующим образом: интеграл в (3.3) разбивается по какой-нибудь формуле численного интегрирования в сумму (для регулярной сетки)

|,J

где Лё(х'га, у'„ ,0) - известное значение поля на сетке, а у - веса формулы интегрирования,

к =_!__

__1_

Ж*,-*:)2-У'.)2+ Полученная система линейных уравнений Ах=у, (у - известные значения Дд(х'т,у'п,0) х- двумерный массив плотности, развернутые в одномерный вектор длины М*Л, где М и N - количества шагов по сетке по х и у соответственно, А - матрица интегрирования размерности МЫ х МЫ) регуляризуется по Лаврентьеву заменой матрицы А на (А+аЕ) и решается любым способом относительно значений ст(х), у^ на сетке.

Опробование на модельных и практических примерах

Заметим, что для определения зон аномальной плотности в горизонтальном слое в Горном институте УрО РАН под руководством В.М. Новоселицкого разработана система «Вектор». При этом определяется некоторый «коэффициент трансформации». Как следует из результатов работ П.С. Мартышко и И.Л. Пруткина, используемая нами методика обладает лучшей разрешимостью при разделении зон аномальной плотности по вертикали. Изложенная методика была опробована на серии модельных примеров. В одном из них по полю, вычисленному от квадратной призмы (ребро в

основании 4000м), отсекаемой двумя параболоидами вращения («купол»), восстанавливалась ее плотность.

Затем задача была усложнена: по полю, вычисленному от квадратной призмы (ребро в основании 4000м), отсекаемой двумя параболоидами вращения («купол»), восстанавливалась ее плотность. Поле было «зашумлено» полем 100 точечных источников, залегавших на глубинах от 200 до 500 метров, максимум суммарного поля которых превосходил на 10% максимум исходного поля (т.е. более 100% «шума»). Затем были последовательно применены описанные в настоящей главе процедуры: исключения боковых источников поля, нахождения поля от источников, залегающих в слое между максимумом верхней границы и минимумом нижней, и решена обратная задача о поиске распределения плотности в слое с «криволинейными» границами. Как видно из рисунков 3.1 и 3.2, аномалия плотности локализована достаточно точно, что позволяет сделать заключение об эффективности применяемой методики.

Рис. 3.1. Исходное (зашумленное) поле Рис. 3.2. Найденное

распределение плотности

С целью дальнейшей проверки и испытания методики была проведена повторная обработка данных по Шершневской и Таманской аномалиям, ранее выполненным в лаборатории математической геофизики ИГФ УрО РАН

(результаты совпали с предыдущей интерпретацией), а также интерпретация данных, предоставленных ЗАО «Гравиразведка» г. Москва (результаты вошли в производственный отчет).

Методики поиска латерального распределения плотности в слоях совместно с Институтом математики и механики УрО РАН (В.ВВасин, Е.Н.Акимова) были реализованы в программах для параллельной вычислительной системы МВС-1000 и успешно опробованы на модельных и практических примерах. Достигнуто значительное ускорение вычислений, что открывает путь к интерпретации значительных массивов данных.

Глава 4. Двумерная задача магниторазведки Постановка задачи.

Рассмотрим следующую постановку двумерной задачи магниторазведки. Пусть в комплексной плоскости в области 1т г<0 находится область О, заполненная веществом с постоянной магнитной проницаемостью ц2- Магнитная проницаемость среды также постоянна и равна ¡1;. Пусть граница 90 области Б - регулярная аналитическая кривая Ь, имеющая горизонтальную асимптоту при 11е(г)-»±оо. В нашем случае можно записать ее уравнение в виде г=з-¡Н+фф, где в - действительный параметр (в частности, Яе(г)=8), <РЙ—а->±<а > 0 . Пусть \^(х,у) - магнитный потенциал первичного поля,

У,(х,у) - внешний, а У2(х,у) - внутренний потенциал области Б. Э" и -внешность и внутренность области О соответственно. Тогда потенциалы V] и У2 с учетом размагничивающего эффекта связаны соотношением

АУк = 0,к = 1,2 -выполняется в Э' и й* соответственно,

У^Уг

дУ2 дУ1 , дЖ

(4.1)

-набО,

а при малых значениях магнитной восприимчивости К

АУк = 0, к = 1,2 - выполняется в и В* соответственно

Введем в рассмотрение функции комплексного переменного (Заморев А. А.)

Эти комплексные функции несут всю информацию о соответствующих действительнозначных функциях потенциала У,(х,у) и могут исследоваться вместо них. Видно, что из (4.1) следует аналитичность и](г) в Б" , и2(г) и и3(г) в Действительно, если У(х,у) - гармоническая функция своих переменных в

д2М д2\ .

некоторой области, то в этой области выполняется—у Л--г- = и.

дх ду

Рассмотрим функцию комплексного переменного и(г), заданную формулой

У1 =у2

(4.2).

Легко видеть, что

аЯе(и(г)) дх

д2У аЬп(и(г))

-7ГТ~п-Г-> и,

ду2 дх

= п

, и, следовательно,

ЭКе(и(2)) _ д1т(и(г))

. Условие

дх ду

ду дх

д2\ д2У

дх

следует из

-=-, и по

дхду дудх

и по теореме о

критерии аналитической функции и(г) является аналитической.

Прямой задачей будем называть задачу поиска функции и^г) при известных О, Ь, к, и3(г), обратной - поиск Б (в частности, границы области) при известных и,(г), и3(г).

Случай без учета размагничивающего эффекта.

Теория двумерной задачи магниторазведки без учета размагничивающего эффекта в случае однородного внешнего поля была развита A.B. Цирульским и Н.В. Федоровой.

Выведем для произвольной гармонической функции V соотношение для

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

да

случае ограниченных объектов подобное соотношение было получено A.A. Чудиновой. Нормаль п к кривой L z=z(s) задается равенством

• z'(s) SV , ,_л

n = i1—^.Поскольку — = (n,gradV), |z(s)| дп

то

— = nxV'+nyV'v = Re дп 1 >

l №1.

Re(-^M)+Im

. V4)

Im(;ru)=

= ;rRe

~Ul

7t UZ -UZ

\z'(s)\) 2 i

Подставляя полученное выражение в (4.1) и пользуясь (4.3), получаем: uxz' -uxz' ~{u2z' -u2z')~ -4kk{u3z' -u3z') (4.4)

Пользуясь свойством

J. z'(s) ,, n ^ '

= Im

показанным A.A. Чудиновой, выражение (4.4) приводим к виду z'(u2 —u,)= -2тпс(и3г' - u3z') (4.5)

Применим интегральную формулу Коши. Интегрирование производится по вещественной оси, и2, и3 не имеют особенностей ниже Г, появляются значения

и3 и и3 на бесконечности по формуле Коши для бесконечных контуров:

— fz'(T)ui(z(T)) _ _2тпс

27С1J T-t 2

r \

•u3(z(t))z'(t) + -

.(4.6)

«у

Т.е. формулой (4.6) задается интегральное уравнение обратной задачи магнитных контактов для малых к.

Для получения решения прямой задачи магнитных контактов воспользуемся так называемым К-уравнением кривой Ь в форме Цирульского. Это

представление в виде 2 = f(z). Видно, что если г=г(8), то

г'О^ГОфуфи Г(г(8)) = ^=1 1ф(8) в нашем случае.

г(б) 1-н<р'(8)

Подставляя в (4.5), получим

щ - и, = -271к(и3 - изг) и применяя формулу Коши, имеем (полагая и3(оо)=0):

2як rf (Qu3(Q , u,(z) = — I dx (4.7).

2m l <;-z

Формулой (4.7) задается решение прямой задачи магнитных контактов. Случай с учетом размагничивающего эффекта

Для ограниченных объектов этот случай был рассмотрен A.B. Цирульским и П.С. Мартшышко.

Для учета размагничивающего эффекта в случае структурных границ перепишем третье равенство из (4.1) в виде

3V,

av,m-n2faw | av,

ön дп ц2 ^ ön Sn

(4.Г)

дЧ

подставляя в (4.Г) выражение для - и пользуясь (4.3), получим после

5п

упрощения

х'(\12 - и!) = - ^^ (щг' -щг' + и3г' - и3г') (4.8) Рассматривая параметр X = ^ , и пользуясь К-уравнением, получим и, =-и, +—-2-1(4.9)

2 Х + 1 1 Х + 1

- прямая задача для магнитных контактов. Из (4.8) получим (подставляя А,):

- (X + 1)г'(з)и2 + Хг'фи, = -г'фи, + ¿фщ - г'(э)и3 (4.10) Интегрируя по Г, получим интефальное уравнение обратной задачи

магнитных контактов

- Ц

'/Ч\ 1 1«

A-f'(s)u,(z(s))ds = U3L __L 2ni * s-t 2 27iif! s-t

, lmt>0. (4.11)

При u3=B (постоянное поле) из (4.11) получаем:

A. fz(s)ui(z(s))¿s_q —L_ fZXs)u1(z(s))ds_5i^>imt>0 2m * s-t 2ni * s-t

При U3 =

z-p

2ni

X_ fZ'(s)ul(z(s))ds = _ J_ rz4s)ul(z(S))ds_ z(t))z,(t) ^^ Inif s-t 2ni f s-t

Вопросы разрешимости. Пусть

+ ¿Л =0 (4.12).

Вычислим явно интегралы в (4.6), это можно сделать, пользуясь теоремой о

f'ix)

вычислении интегралов типа Коши для функций вида J—^-у (p(x)dT:

dx =

_ Yk .

27ri/(z(t)-akXx-t)"b t-tk'

Г _

где г^^аь к=1,...,М и ^с^Ьь к=1,...,п;

В случае без учета размагничивания, когда у функции и3 существует первообразная Г, можно записать дифференциальное уравнение

__Цд

2-t_r + Ыч —wi

dt = -2racu3 dz, (4.13)

откуда

k=l k=l j-1 j\t-ckJ

= -27tKF(z),

z = t-ih—— F"1 2 пк

ЕпшЮ-ЁЕ^Цг

(4.14).

в случае с учетом размагничивания уравнение примет вид

2як

{tnlnb^-i^y^rl^I^Hni'-O-ii^y)

*.i j\f-cj l.»-' Jv - ct У )t

ih + t

(4.15)

Таким образом, доказана

Теорема. Если и^г) представляется в виде (4.12), и3(г) имеет первообразную Р, то решение обратной задачи представляется без учета размагничивания в виде (4.14), или с учетом размагничивания в виде (4.15).

Решение обратной задачи сводится к решению нелинейной, вообще говоря, системы (4.14) или (4.15) при подстановке в нее соотношений 2^к)=ак, к=1,...,Ы и г(ск)=Ьк, к=1,...,п относительно ^ Ск.

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

1. Поле и^г) подбирается полем системы стержней, имеющих поле вида

либо в виде (4.12), при этом решается задача минимизации нелинейного функционала.

2. По формуле (4.14) решается задача о нахождении гф для случая магнитных контактов, при этом потребуется решать нелинейные уравнения для нахождения ^ ск.

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

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

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

осуществить решение задачи поиска распределения избыточной плотности в горизонтальном слое и в слое с «криволинейными» границами;

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

осуществить двух- и трехмерную визуализацию исходного поля и изучаемых структур. Программа написана на языке Object Pascal в среде программирования Delphi 5 с использованием API OpenGL для осуществления визуализации. Программа откомпилирована под ОС MS Windows ХР (тестировалась также под MS Windows 98).

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

Программа состоит из двух трех блоков - пересчета поля, решения обратной задачи и визуализации. Основным алгоритмом блока пересчета поля является алгоритм, изложенный в главе 2. Для пересчета поля вверх использовался метод интегрирования по составной формуле трапеций, а для пересчета поля вниз использовались 2 метода обращения (с регуляризацией) матрицы

интегрирования: точный (Гаусса) и приближенный (метод наискорейшего спуска). Как показали эксперименты, при пересчете поля вниз (с регуляризацией) невязка в 1% достигается менее, чем за 60 итераций, что занимает 18 секунд, а непосредственное обращение матрицы занимает 140 секунд (Процессор Intel Celeron 2.6 ГГц). Для решения двумерной задачи Дирихле использовался метод сеток с обращением матрицы методом Гаусса. Решение требовало менее 8 секунд.

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

Рис. 5.1. Основное окно программы. Блок визуализации осуществляет графическое представление информации о поле в двух видах: в виде трехмерной поверхности г=А[х,у) (визуализация

средствами OpenGL) и в виде двухмерной градиентной картинки (оттенки серого соответствуют значениям поля).

Для градиентной картинки предусмотрена возможность построения системы изолиний (методом разделения на указанное число уровней яркости). Имеется возможность строить разрезы (двумерные графики значений функции) вдоль профилей, указываемых пользователем.

Программа использовалась для интерпретации данных гравиразведки по 3 регионам.

ЗАКЛЮЧЕНИЕ И ОСНОВНЫЕ ВЫВОДЫ

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

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

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

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

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

5. Алгоритмы решения обратной задачи о распределении плотности в плоском и криволинейном слое реализованы в компьютерной программе и были применены при обработке реальных геофизических измерений (произведена обработка гравитационных данных по 3 аномалиям). В соавторстве с П.С. Мартышко результаты были представлены в 8 публикациях.

6. Совместно с Институтом математики и механики УрО РАН (Е.Н Акимова, В.В. Васин) осуществлено (для модельных и практических примеров) опробование решения задачи о распределении плотности в криволинейном слое на многопроцессорной параллельной вычислительной системе МВС-1000 (язык программирования - MPI Fortran), доказана эффективность

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

Сведения о публикациях.

Результаты исследований по теме диссертации опубликованы в следующих работах:

1. Мартышко П.С., Пруткин ИЛ., Кокшаров Д.Е. Об одном алгоритме построения плотностных разрезов по гравитационному полю. // Доклады VI междисциплинарного международного симпозиума «Закономерности строения и эволюции геосфер». Хабаровск, 2004.

2. Мартышко П.С., Кокшаров Д.Е. О построении плотностных разрезов по гравитационному полю. // III Всероссийская конференция «Алгоритмический анализ неустойчивых задач». Екатеринбург, 2004.

3. Martyshko P.S., Koksharov D.E.. On the construction of density sections using gravity data. // 66th EAGE Conference & Exhibition. Extended abstracts. Paris, 2004.

РОС. НАЦИОНАЛ ! „ БИБЛИОТЕКА С Петербург W IM шаг

4. Мартышко П.С., Кокшаров Д.Е. Двумерные задачи магниторазведки для магнитных контактов. // Проблемы теоретической и прикладной математики: Труды 36-й Региональной молодежной конференции. Екатеринбург: УрО РАН, 2005.

5. Мартышко П.С., Кокшаров Д.Е. Об одном методе решения обратных задач гравиметрии. // Глубинное строение, геодинамика, мониторинг, тепловое поле Земли, Интерпретация геофизических полей. Третьи научные чтения Ю.П. Булашевича. Материалы конференции. ИГФ УрО РАН, Екатеринбург, 2005.

6. Martyshko P.S., Koksharov D.E. New approach to the density model construction using gravity data. // Regularities of the Structure and Evolution of Geospheres: Materials of VII International Interdisciplinary Scientific Symposium. Vladivostok, 20-24 September, 2005.

7. Кокшаров Д.Е. Компьютерная технология обработки данных гравиразведки для решения задачи о распределении избыточной плотности в слое. // Геофизика-2005. Пятая международная научно-практическая геолого-геофгоическая конференция-конкурс молодых ученых и специалистов. Тезисы докладов. - СПб.: СПбГУ, ВВМ, 2005.

8. Мартышко П.С., Кокшаров Д.Е. Об определении плотности в слоистой среде по гравитационным данным. // Геофизический журнал, Т. 27. №4, 2005.

i í

»24201

РНБ Русский фонд

2006-4 26768

I1'

Содержание диссертации, кандидата физико-математических наук, Кокшаров, Дмитрий Евгеньевич

ВВЕДЕНИЕ.

1. МЕТОДЫ РЕШЕНИЯ ОБРАТНЫХ ЗАДАЧ ГРАВИМЕТРИИ.

2. ВЫДЕЛЕНИЕ АНОМАЛИЙ ПОЛЯ В УКАЗАННОЙ ОБЛАСТИ,

4 ПОДГОТОВКА ДАННЫХ.

2.1. Выделение аномалий поля в указанной области.

2.2. Подготовка данных.

3. ОБРАТНАЯ ЗАДАЧА ГРАВИМЕТРИИ.

3.1. Структурные границы. Метод локальных поправок.

3.2 Слой переменной плотности.

3.3 Опробование на модельных и практических примерах.

4. ДВУМЕРНАЯ ЗАДАЧА МАГНИТОРАЗВЕДКИ.

4.1. Постановка задачи.

4.2. Случай без учета размагничивающего эффекта.

4.3. Случай с учетом размагничивающего эффекта. ф 4.4. Вопросы разрешимости.

4.5. Решение прямой задачи для предельных случаев.

5. КОМПЬЮТЕРНАЯ ПРОГРАММА ОБРАБОТКИ ДАННЫХ ГРАВИРАЗВЕДКИ.

5.1. Общее описание программы.

5.2. Описание интерфейса программы.

Введение Диссертация по наукам о земле, на тему "Алгоритмы и новые компьютерные технологии решения структурных обратных задач гравиметрии и магнитометрии"

Обоснование актуальности темы

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

В теории потенциальных геофизических полей различают два важнейших класса задач - прямые и обратные. Прямая задача - нахождение поля, порождаемого геологической структурой с известными свойствами (параметрами) - границами, распределением плотности или намагниченности. Такие задачи (моделирования полей) достаточно хорошо исследованы как для гравитационных, так и для магнитных полей в двух- и трехмерном случае (можно отметить работы 4 В.Н. Страхова, Ю.И. Блоха, П.С. Мартышко, С.М. Оганесяна, В.И. Старостенко, а также М. Talwani, S.P. Sharma, Р. Weidelt и других). Обратные задачи -определение физических параметров пород по измеренным полям. При разработке теории решения обратных задач гравиметрии и магнитометрии лучше всего исследован двумерный случай. Результаты В.Н. Страхова, A.B. Цирульского, В.К. Иванова, П.С. Мартышко, Н.В. Федоровой, Ф.И. Никоновой позволяют эффективно находить решение задачи -эквивалентное семейство решений (существенный прогресс в этом направлении достигнут с применением результатов ТФКП). В трехмерном случае обратные задачи решаются менее успешно. Причин этому несколько. Во-первых, такие задачи сводятся к решению обратных задач математической физики, а последние, в свою очередь, являются некорректно поставленными, что существенно затрудняет интерпретацию полей. Общая теория и методы решения некорректно Ь поставленных задач были разработаны А.Н. Тихоновым, В.К. Ивановым, М.М. Лаврентьевым, а применительно к задачам геофизики методы регуляризации разрабатывали В.Н. Страхов, В.Б. Гласко. Следует отметить работы А.И. Кобрунова, посвященные теоретическим основам решения обратных задач геофизики в трехмерных постановках и В.М. Новоселицкого, много сделавшего для практических реализаций методик интерпретации; интересные практические методы разработаны В.Н. Страховым, В.И. Старостенко, П.С. Мартышко, И.Л. Пруткиным, A.C. Долгалем. Во-вторых (это следует из результатов по теории эквивалентности В.Н. Страхова, A.B. Цирульского, С.М. Оганесяна, В.И. Старостенко и других), решить однозначно обратную задачу при интерпретации геофизических полей естественной природы в принципе невозможно. Кроме того, задачи такого класса требуют, как правило, существенных вычислительных ресурсов, что до последнего времени было недоступно. Следовательно, решение задач практической интерпретации данных потенциальных геофизических полей остаются актуальными, особенно с привлечением современных вычислительных средств.

Связь работы с научными программами, планами, темами

Работа выполнялась в соответствии с плановой тематикой НИР Института по теме 01.200.2 09053.

Цель и задачи исследования.

Предмет и объект исследования.

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

Научная новизна полученных результатов.

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

В рамках одного пакета осуществляется полный цикл обработки данных:

• загрузка исходных файлов;

• исключение боковых (находящихся вне рассматриваемой области) источников поля;

• выделение поля от аномалий, залегающих в изучаемом слое;

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

3. Совместно с ИММ УрО РАН (В.В. Васин, E.H. Акимова) предложенные алгоритмы решения структурной трехмерной обратной задачи (в линейной постановке) были реализованы в технологии интерпретации гравиметрических данных с использованием параллельных вычислений, реализована и опробована на модельных примерах и практических данных программа для кластера МВС-1000.

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

Практическая значимость полученных результатов.

Компьютерная методика решения линейной трехмерной обратной задачи для плоских и слоев, границы которых являются поверхностями, имеющими плоскую горизонтальную асимптоту, была успешно применена для интерпретации практических данных по 3 регионам, результаты были переданы заказчикам и включены в отчеты ЗАО «Пермь-Лукойл» (г. Пермь) и ЗАО «Гравиразведка» (г. Москва).

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

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

Сведения об апробации результатов диссертации. Результаты научной работы доложены:

1. 66th EAGE Conference & Exhibition. Extended abstracts. Paris, 2004.

2. VI междисциплинарного международного симпозиума. Хабаровск, 2003.

3. 36-й Региональная молодежная конференция «Проблемы теоретической и прикладной математики», ИММ УрО РАН, 2005.

4. Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей. Международный Семинар им Д.Г. Успенского, 2004.

5. 6 Уральская молодежная научная школа по геофизике. Горный институт УрО РАН. Пермь, 2005.

Ю.П. Булашевича. 2005.

Сведения о публикациях.

Результаты исследований по теме диссертации опубликованы в следующих работах:

1. Мартышко П.С., Пруткин И.Л., Кокшаров Д.Е. Об одном алгоритме построения плотностных разрезов по гравитационному полю. // Доклады VI междисциплинарного международного симпозиума «Закономерности строения и эволюции геосфер». Хабаровск, 2004.

3. Martyshko P.S., Koksharov D.E. On the construction of density sections using gravity data. // 66th EAGE Conference & Exhibition. Extended abstracts. Paris, 2004.

7. Кокшаров Д.Е. Компьютерная технология обработки данных гравиразведки для решения задачи о распределении избыточной плотности в слое. // Геофизика-2005. Пятая международная научно-практическая геолого-геофизическая конференция-конкурс молодых ученых и специалистов. Тезисы докладов. - СПб.: СПбГУ, ВВМ, 2005.

8. Мартышко П.С., Кокшаров Д.Е. Об определении плотности в слоистой среде по гравитационным данным. // Геофизический журнал, №4. Институт Геофизики HAH Украины. Киев, 2005.

Работа частично (компьютерная технология в методе поиска структурных границ) поддержана грантом УрО РАН для молодых ученых за 2005 г.

Автор выражает глубокую благодарность научному руководителю - д.ф.-м.н., профессору П.С. Мартышко, а также д.ф.-м.н. И.Л. Пруткину, чл.-корр. РАН, профессору В.В. Васину, к.ф.-м.н. E.H. Акимовой, к.ф.-м.н. Н.В. Федоровой, к.ф.-м.н. Ф.И. Никоновой за консультации и интересные обсуждения в процессе работы.

Заключение Диссертация по теме "Геофизика, геофизические методы поисков полезных ископаемых", Кокшаров, Дмитрий Евгеньевич

Заключение

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

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

Библиография Диссертация по наукам о земле, кандидата физико-математических наук, Кокшаров, Дмитрий Евгеньевич, Екатеринбург

1. Андреев В. И. Моделирование геологических образований методами пространственной гравиметрии. — М: Недра, — 1992. —- 224с.

2. Байков В.А., Жибер A.B. Уравнения математической физики. Москва-Ижевск: Институт компьютерных исследований, 2003, 256 с.

3. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Лаборатория базовых Знаний, 2002. 632 с.

4. Булах Е. Г Автоматизированная система интерпретации гравитационных аномалий. — К.: Наук, думка, 1973. — 111с.

5. Булах Е.Г., Ржаницин В.А., Маркова М.Н. Применение метода минимизации для решения задач структурной геологии по данным гравиразведки. К.: Наук, думка, 1976. - 220 с.

6. Буллах Е.Г., Маркова М.Н. Метод простого моделирования при решении задач гравиметрии в классе трехмерных контактных поверхностей. // Геофизический журнал. 2003. - №1, Т. 25, С. 92-100.

7. Великович А. Л., Зельдович Я. Б. Об одном подходе к решению обратной задачи теории потенциала.- Докл. АН СССР, т. 212, № 3, с. 580-583.

8. Вере Л., Джон Ф., Шехтер М. Уравнения с частными пропзводпымп. М.: Мир. 1966. 352 с.

9. Гахов Ф.Д. Краевые задачи. М., Физматгиз, 1963. 640 е., ил.

10. Голов И.Н, Сизиков B.C. Обратная задача гравиметрии с использованием сфероидов, нелинейного программирования и регуляризации. // Геофизический журнал, №3, 2005.

11. Гольцман Ф.М., Калинина Т.Б. Статистическая интерпретация магнитных и гравитационных аномалий. -Ленинград: Недра, 1983. 245 с.

12. Иванов В. К. Об определении гармонических моментов возмущающих масс по производной гравитационного потенциала, заданной на плоскости.- Изв. АН СССР. Сер. геогр. и геофиз., 1950, т. 14, № 5, с. 403414.

13. Иванов В.К., Васин В.В., Танана В.П. Теория линейных некорректных задач и ее приложения. М., «Наука», 1978, 206 с.

14. Иоффе А.Д., Тихомиров В.М. Теория экстремальных задач. М.: Наука, 1974. -479 с.

15. Кобрунов А.И., Варфоломеев В. А. К построению пространственных распределений масс с нулевым гравитационным эффектом // Геофизический журнал. 1983. - Т.5. - N 3. С. 53-56.

16. Кобрунов А. И. Теоретические основы решения обратных задач геофизики. — Ухта: Изд-во. УИИ, 1995. — 228 с.

17. Кобрунов А. И., Варфоломеев В А. Об одном методе е-эквпвалентиых перераспределений и его использовании при интерпретации гравитационных полей.-Изв. АН СССР. Физика Земли, 1981, № 10, с. 2544.

18. Кобрунов А.И. К вопросу об интерпретации аномальных гравитационных полей методом оптимизации (трехмерная задача) // Изв. АН СССР.Физика Земли. 1979.- N 10. - С. 67-77.

19. Кобрунов А.И. О классах оптимальности решений обратнойзадачи гравиразведки // Изв.АН СССР. Физика Земли. 1982. -N2.-0. 100-107.

20. Кобрунов А.И. О построении решений обратной задачи гравиразведки в классе распределений плотности // Докл. АН УССР,Б. 1977. - N 12. - С. 1077-1079.

21. Кобрунов А.И. Разрешимость и эквивалентность в обратной задаче гравиразведки для нескольких плотностных границ // Изв. АН СССР. Физика Земли. 1983.- N 5.- С. 67-75.

22. Колмогоров А.Н., Фомин С. В. Элементы теории функций и функционального анализа. М.: Наука, 1976. - 542 с.

23. Корнейчук Н.П. Экстремальные задачи теории приближения. М.: Наука, 1976.-320 с.

24. Красовский С. С. Отражение динамики земной коры континентального типа в гравитационном поле. — Киев: Наук, думка, 1981. — 264с.

25. Лаврентьев М.М. и др. Некорректные задачи математической физики и анализа. М.: Наука, 1980. - 286 с.

26. Лаврентьев М.М. О некоторых некорректных задачах математической физики. Новосибирск, 1962.

27. Ладыженская О. А. Краевые задачи математической физики. М.: Наука, 1973,408 с.31 .Маловичко А. К. Об интерпретации гравиметрических наблюдений в связи с поисками структур, перспективных на нефть и газ // Прикл. геофизика. — 1955. — Вып. 13. — С. 36—79.

28. Маргулис А. С. Гармонические плотности и обратные задачи потенциала-В кн.: Теория и методика интерпретации гравимагнитных полей. Киев: Наук, думка, 1981, с. 130-136.

29. Маргулис A.C. Гармонические плотности и обратные задачипотенциала// Теория и методика интерпретации гравимагнитных полей. Киев: Наукова ' думка, 1981.-С. 130-136.

30. Маргулис A.C. К теории потенциала в классах ЬР(П)//Изв. вузов. Математика. 1982. - N 1. - С. 33-41.

31. Маргулис A.C. Теория потенциала для плотностей класса Ьр(П)и ее применение к обратным задачам гравиметрии // Теория и практика интерпретации гравитационных и магнитных аномалий. -Киев: Наукова думка, 1983.-С. 188-197.

32. Мартышко П.С., Кокшаров Д.Е. Двумерные задачи магниторазведки для магнитных контактов. // Проблемы теоретической и прикладной математики: Труды 36-й Региональной молодежной конференции. Екатеринбург: УрО РАН, 2005.

33. Мартышко П.С., Кокшаров Д.Е. О построении плотностных разрезов по гравитационному полю. // Тезисы докладов всероссийской конференции «Алгоритмический анализ неустойчивых задач». Екатеринбург, 2004.

34. Мартышко П.С., Пруткин И.Л. Технология разделения источников гравитационного поля по глубине. // Геофизический журнал. Т.25. № 3, 2003, с. 159-168.

35. Мартышко П.С., Пруткин И.Л., Кокшаров Д.Е. Об одном алгоритме построения плотностных разрезов по гравитационному полю. «Строение и динамика геосфер», Хабаровск, 2004.

36. Мартышко П.С., Пруткин И.Л., Начапкин Н.И., Шестаков А.Ф. Отчет по теме «Разработка технологии разделения источников гравитационного поля по глубине». ИГФ УрО РАН. Екатеринбург, 2000.

37. Мартышко П.С., Пруткин И.Л., Шестаков А.Ф. Отчет по теме «Вопросы математической интерпретации гравитационных данных по Соликамской впадине». ИГФ УрО РАН. Екатеринбург, 2001.

38. Мартышко П.С., Пруткин И.Л., Шестаков А.Ф. Отчет по теме «Вопросы математической интерпретации гравитационных данных. II. (по Соликамской впадине)» ИГФ УрО РАН. Екатеринбург, 2002.

39. Михайлов В. П. Дифференциальные уравнения в частных производных. М.: Наука, 1983. 424 с.

40. Новоселицкий В.М., Маргулис A.C. Фурье- аналогия в обратных задачах гравиметрии для некоторых "нефтяных" и "планетарных" плотностных моделей // Изв. АН СССР. Физика Земли. 1982. - №4. - С. 68 - 82.

41. Оганесян С. М., Старостенко В.И. L- псевдорешения и их использование для построения тел с нулевым внешним гравитационным полем// Изв. АН СССР.Физика Земли. 1984. - N 2. - С.51-62.

42. Оганесян С. М., Старостенко В.И. Некоторые свойства L -псевдорешений и элемента, минимизирующего обобщенный функционалА.Н. Тихонова// Докл. АН УССР. 1982. - т. 263. - N 3. - С.540-542.

43. Оганесян С.М., Старостенко В.И. Тела нулевого внешнего гравитационного потенциала: о забытых работах и современном состоянии теории. // Физика Земли. 1985. №3. С. 49-62.

44. Оганесян СМ. Решение обратной задачи гравиметрии в классераспределений плотности // Докл. Ан УССР, Б. -1981. N 6. -С.39-43.

45. Пруткин И.Л. О предварительной обработке измерений, заданных на площади // Методы интерпретации и моделирования геофизических полей. Свердловск: УрО АН СССР, 1988. С. 11-15.

46. Соболев С. Л. Введепие в теорию кубатурпых формул. М.: Наука, 1974. 808 с.

47. Старостенко В. И., Черная Н. Н., Черный А. В. Обратная задача для контактных поверхностей 1,11,III // Физика Земли. — 1992. — №6. — С. 48— 58; —1993. — №7. — С. 47—66.

48. Старостенко В.И. Устойчивые численные методы в задачах гравиметрии. -Киев: Наук, думка, 1978. 228 с.

49. Старостенко В.И., Исаев В.И., Пятаков Ю.В. Решение обратной задачи для контактов осадочных пород. // Геофиз. журн. 1993. - №1. - С. 62-71.

50. Страхов В. Н. К теории метода подбора // Изв. АН СССР, сер. геофиз. — 1964. — №4. — С. 494— 509.

51. Страхов В.Н. О новом этапе в развитии теории интерпретации гравитационных и магнитных аномалий // Изв. АН СССР. Физика Земли. -1977.-К2-С.20-41.

52. Страхов В.Н. Об общих решениях обратных задач гравиметрии и магнитометрии // Изв. вузов. Геология и разведка. 1978. - N 4. -С. 104-117.

53. Страхов В.Н. Эквивалентность в плоской задаче гравиметрии при переменной плотности масс// Изв. АН СССР.Физика Земли. -1977. N 5. - С. 48-60.

54. Страхов В.Н., Гольшмидт В.И., Калинина Т.Б., Старостенко В.И. Состояние и перспективы развития в СССР теории интерпретации гравитационных имагнитных аномалий // Изв. АН СССР. Физика Земли,- 1982.- N 5. С. 1133.

55. В.К. Иванов и др. Теория линейных некорректных задач и ее приложения. -М.: Наука, 1978.-206 с.

56. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1979. 285 с.

57. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. М: Наука. Гл. ред. физ.-мат. лит., 1990. 232 с.

58. Тихонов А.Н., Леонов А.С., Ягола А.Г. Нелинейные некорректные задачи. -М.: Наука. Физматлит, 1995. 312 с.

59. Торге В. Гравиметрия: Пер. с англ. М., Мир. 1999. - 429 е., ил.

60. Хилле Э., Филлипс Р. Функциональный анализ и полугруппы. -М.: Изд. ин. лит., 1962. 829 с.

61. Цирульский А.В. Функции комплексного переменного в теории и методах потенциальных геофизических полей. Свердловск: УрО АН СССР, 1990.

62. Цирульский А.В., Мартышко П.С. Об учете размагничивающего эффекта в задачах магниторазведки // Изв. АН СССР. Физика Земли. 1979. №3.

63. Чудинова А.А. Интегральное уравнение обратной задачи электроразведки. // Изв. вузов. Сер. матем. 1965. №6. С. 1080-1088.

64. Шалаев С. В. Геологическое истолкование геофизических аномалий с помощью линейного программирования. Л.: Недра, 1972. - 144 с.

65. Akima, Н. 1978. A method of bivariate interpolation and smooth surface fitting for irregularly distributed data points. ACM Trans, on Math. Software 4: 148159.

66. Barbosa V. C. F., Silva J. В. C., Medeiros W. E., 1997, Gravity inversion of basement relief using approximate equality constraints on depths: Geophysics, 62, 1745-1757.

67. Botl M H. P., 1960, The use of rapid digital computing methods for gravity interpretation of a sedimentary basin: Geophysical Journal of the Royal Astronomical Society, 3, 63-67.

68. Briggs I.C. 1974. Machine contouring using minimum curvature. Geophysics 39: 39-48.

69. Cooke R.A. Mostaghimi S. 1992. A microcomputer based routine for obtaining mean watershed precipitation from point values. Computers and Geosciences 18: 823-837.

70. Cordell L., Henderson R. G., 1968, Iterative three dimensional solution of gravity anomaly data using a digital computer: Geophysics, 33, 596-601.

71. Davis J. C. 1986. Statistics and data analysis in geology (2nd ed.). New York: John Wiley & Sons.

72. Franke R, Nielson G. Smooth Interpolation of Large Sets of Scattered Data . InternationalJournal for Numerical Methods in Engineering, 1980.

73. Guillen A., Menichelti V., 1984, Gravity and magnetic inversion with minimization of a specific functional: Geophysics, 49, 1354-1360.82.1ngber L., 1989, Very fast simulated re-annealing: Journal of Mathematical and Computer Modelling, 12, 967-973.

74. Kelway P. S. 1974. A scheme for assessing the reliability of interpolated rainfall estimates. Jour. Hydrology 21: 247-267.

75. Kirkpartrick S. G., Delatt C. D. Jr., Vecchi M. P., 1983, Optimization by simulated annealing: Science, 220, 671 -680.

76. Lam N. S. 1983. Spatial interpolation methods review. The American Cartographer 10: 129-149.

77. Last B. J., and K. Kubik, 1983, Compact gravity inversion: Geophysics, 46, 713721.

78. Lawson C.L. 1972. Generation of a triangular grid with application to contour plotting. Tech. Memo. 299, Sect. 914, Jet Propulsion Lab., Caltech, Pasadena, California.

79. Leao J. W. D., P. T. L. Menezes, J. F. Beltrao. and J. B. C. Silva, 1996, Gravity inversion of basement relief constrained by knowledge of depth at isolated points: Geophysics, 61, 1702-1714.

80. Martyshko P.S., Koksharov D.E. On the construction of density sections using gravity data. Extended Abstracts. 66th EAGE Conference & Exhibition. Paris, 2004.

81. Mauriello P., Patella D. Localization of maximum-depth gravity anomaly sources by a distribution of equivalent point masses. Geophysics. Vol. 66, No. 5, 2001. P.1431-1437.

82. McBratney A.B., Webster R. 1986. Choosing functions for semi-variograms of soil properties and fitting them to estimates. Jour. Soil Science 37: 617-639.

83. Montefusco L.B. Casciola G. 1989. CI surface interpolation. ACM Transactions on Mathematical Software 15: 365-374.

84. Nagihara, S., Stuart A. H. 2001. Three dimensional gravity inversion using simulated annealing: Constraints on the di-apiric roots of allochthonous salt structures: Geophysics, 66, 1438-1449.

85. Preusser A. 1990. CI- and C2- interpolation on triangles with quintic and nonic bivariate polynomials. ACM Transactions on Mathematical Software 16: 253257.

86. Sen M. K., and P. L. Stoffa, 1995. Global optimization methods in geophysical inversion: Elsevier Science Publ.

87. Sen M. K., and P. L. Stoffa, 1996. Bayesian inference, Gibbs' sampler and uncertainty estimation in geophysical inversion: Geophysical Prospecting, 44, 313-350.

88. Talwani M., Ewing M. Rapid computation of gravitational attraction of three-dimensional bodies of arbitrary shape. Geophysics, vol. XXV, No 1, 1960. pp. 203-255.

89. Talwani, M., Worzel J. L., Landisman M. 1959, Rapid gravity computations for two-dimensional bodies with application to the Mendocino submarine fracture zone: Journal of Geophysical Research, 64, 49-59.

90. Tarantola, A., 1987, Inverse problem theory: Methods of data fitting and model parameter estimation: Elsevier Science Publ.

91. Tevis J.W., Whittaker A.D., McCauley D.J. 1991. Efficient use of data in the kriging of soil pH. ASAE paper 91-7074. St. Joseph, MI: ASAE.

92. Watson D.F., Philip G.M. 1985. A refinement of inverse distance weighted interpolation. Geo-Processing 2: 315- 327.

Информация о работе

Кокшаров, Дмитрий Евгеньевич
кандидата физико-математических наук
Екатеринбург, 2005
ВАК 25.00.10

Диссертация

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

Автореферат

Похожие работы