ФИЗИКА ЗЕМЛИ, 2008, № 7, с. 21-27
УДК 550.8
ОБРАТНЫЕ ЗАДАЧИ ГРАВИМЕТРИИ ДЛЯ СОВОКУПНОСТИ ТЕЛ КЛАССА Л.Н. СРЕТЕНСКОГО
© 2008 г. Е. Г. Булах, М. Н. Маркова
Институт геофизики НАН Украины им. С И. Субботина, г. Киев Поступила в редакцию 15.01.2007 г.
Решается обратная задача гравиметрии, когда искомая геологическая модель есть совокупность тел класса Л.Н. Сретенского. Для каждого тела фиксировано положение средней плоскости. Нужно найти положения кровли и подошвы тела. Эти поверхности аппроксимируются аналитическими функциями и определяются параметрами. Решение задачи иллюстрируется примером.
PACS: 91.10.0p
1. ВВЕДЕНИЕ
Обратные задачи гравиметрии и магнитометрии чаще всего решаются в фиксированных классах геологических моделей. В рудной геофизике приходится иметь дело не только с одиночными телами. Аномальное поле обычно обусловлено совокупностью тел. Это обстоятельство должно быть учтено при разработке алгоритмических и программных решений обратных задач.
Аппроксимационный подход становится ведущим при решении интерпретационных задач гравиметрии и магнитометрии. Развитие этого подхода дано в большом числе работ В.Н. Страхова (1973, 1978, 2005 и др.).
Л.Н. Сретенский (1954) установил один класс геологических моделей. Введено понятие средней плоскости. Будем говорить, что тело обладает средней плоскостью, если существует такая плоскость Р, называемая средней, когда всякая прямая, перпендикулярная к этой плоскости, может пересекать поверхность тела только в двух точках, расположенных по разные стороны от плоскости Р. Теоремой утверждается, что если геологическое тело обладает средней плоскостью, его центр тяжести находится внутри тела, избыточная плотность постоянная, то обратная задача ньютонового потенциала имеет единственное решение.
Таким образом, тело, которое относится к классу Л.Н. Сретенского, обладает средней плоскостью, верхней и нижней граничными поверхностями. Каждая поверхность есть однозначная функция координат точек. Можно по-разному подходить к аппроксимационному построению этих граничных поверхностей. Будут создаваться разные подклассы геологических объектов. В этих подклассах геологическая модель параметризуется. Устанавливается такая совокупность параметров, численные значения которых опре-
деляют местоположение и размеры тела. Два таких аппроксимационных построения были рассмотрены ранее - Булах Е.Г. (2002). Один из них будет реализован в этой работе.
2. АППРОКСИМАЦИЯ МОДЕЛИ
Выберем систему координат. Пусть начало координат расположено в точке дневной поверхности. Ось аппликат ОЪ или 0£, направлена вертикально вниз. Координатная плоскость ХОУ или %Оп - горизонтальная. Она совпадает с дневной поверхностью, если последняя также горизонтальная.
Примем такое правило. Если точка М принадлежит гравитирующему телу, то будем писать М(%, п, С). Если же точка М внешняя, то запишем М(х, у, z). Пусть задано некоторое тело. Оно ограничено замкнутой поверхностью и заполнено массами, избыточная плотность которых постоянная. Зафиксируем горизонтальную плоскость z = Н0. Она пересекает аномальное тело и является для него средней плоскостью. В сечении выделена область Б, принадлежащая гравитирующему объекту. Пусть эта область ограничена прямоугольным контуром. В этом случае постулируем
Б : [а < % < Ь; с < п < ^ (1)
Гравитирующее тело обладает верхней и нижней поверхностями. Запишем
Ъ = Но - Ъ1 (%,п); Ъ = Но + Ъ2(%,п)
mF,
Z1 = X 4
•( 1)
sin
пt(£, - a) . пt(n - c)
b - a
sin-
d - c
t = 1 mF
Z 2 = X 4
,(2)
sin
пt(£, - a) . пt(n - c)
b - a
sin-
d - c
t = 1
Здесь введены две функции. Они однозначны, внутри области D должны быть положительными и на границе этой области обращаются в нуль. Модель будет рассматриваться последовательным увеличением числа слагаемых.
Если геологическая модель состоит из m тел, то она описывается такой совокупностью параметров
P = {m; [а; H„; (a, b, c, d);
(mFt = 1, 2, ...mFi);
(mF2, 2), t = 1, 2, ...mF2)]j, j = 1, 2, ...m }
(3)
3. ПРЯМАЯ ЗАДАЧА ГРАВИМЕТРИИ
Выше была описана гравитирующая модель. Пусть она содержит одно тело. Аномалия силы тяжести, обусловленная избыточными массами тела, запишется так
Аg (x, y, z) = -kа x
b b
x
я
[(5 - x) 2 + (n - y )2 + (Z - z )2 ]1/2
H + Z 2
(4)
d5 dn.
H„- Z1
bb
Аg (x, y, z) = kа
U ö1
x
x
1-1-
(H „ + Z 2- z )2- (H„- Z1- z )2"
Q2
x d5dn
(4.1)
А1. Подготовка к решению обратной задачи.
Обычно, когда решаются задачи рудной геофизики, аномальное поле может быть осложнено сторонними влияниями. Говорят, что поле имеет фоновую составляющую. На первом этапе интерпретационных работ исходное поле может быть аппроксимировано гармонической функцией. При этом фоновое влияние учитывается дополнительным слагаемым.
Заменим поле (5) полем модельных источников - например, совокупностью шаровых тел или совокупностью трехстержневых тел [Булах, 1967; 2000]. Если поле (5) было осложнено фоновым влиянием, то эта фоновая составляющая аппроксимируется несложной линейной или квадра-тической алгебраической функцией. Из исходного поля фоновое влияние исключается. Теперь нетрудно получить целую серию аномальных гравитационных полей. Пусть это будет такой набор полей
Agn(x, y, z) ^ Un(x, y) = = {Agn(x, y), 5Agn(x, y), Vnzz(x, y), V^zz(x, y)}
(6)
Здесь уже выполнено первое интегрирование по аппликате. Под знаком интеграла разность двух функций, которые могут быть достаточно близки между собой. Для уменьшения вычислительной погрешности преобразуем это выражение таким образом
Здесь записано следующее. Поле, заданное в точках на дневном рельефе, преобразуется в поле на горизонтальной плоскости г = 0. Из исходного поля легко получить поле вариации аномалии силы тяжести относительно поля в фиксированной точке. Исходное поле можно продифференцировать и получить поля высших производных.
Описанный прием позволяет повысить результативность качественной интерпретации. Высшие производные потенциала хорошо локализуют объекты. Анализ таких полей позволяет установить число локальных источников и координаты их эпицентров. В интерпретационный процесс может быть включено любое исходное поле.
А2. Решение задачи. Для решения обратной задачи может быть выбрано какое-либо исходное поле, приведенное в записи (6). В этом поле зафиксируем п точек и получим первый массив:
Un(xi, y) = Un(i), i = 1, 2, ... n.
(7)
Q12 = (5-x)2 + (n -y)2 + (H„-Z1 -z)2; Q22 = (5 - x)2 + (n - y)2 + (H„ + Z2- z)2.
4. ОБРАТНАЯ ЗАДАЧА
Задано поле аномалии силы тяжести. В этом поле выбрано n точек и получен такой массив исходных данных
Agn(xi, yi, zi) = Agn(i), i = 1, 2, ...n. (5)
Теперь на основании исходного поля и априорных данных интерпретатором составляется геологическая модель. Она описывается совокупностью параметров - запись (3). Естественно, что численные значения некоторых параметров нужно постулировать заранее. Их значения не будут изменяться. Другие параметры получают начальные значения, и эти величины будут изменяться в итерационном вычислительном цикле. По внешнему аномальному полю можно определить число аномальных тел. Из априорных данных установим положение средних плоскостей. В каждой средней плоскости выделим
a c
a c
У
X
Рис. 1. Поле аномалии силы тяжести в точках дневной поверхности.
область D - запись (1). Таким образом, запись (3) преобразуется так
(8)
р = [ w 1; W2 ];
№ 1 = [т; (Но, а, Ь, с, й), 1 = 1, 2, ...т]; №2 = [а; (т^1; S<t1), г =1, 2, ...т^1);
(т^2; £(2), г = 1, 2, ...тр)1 = 1, 2, ...т].
Здесь последовательность Ш объединяет постоянные параметры, а последовательность №2 - переменные. Все параметры должны получить начальные численные значения. Может быть решена прямая задача и получено теоретическое поле
иг(х, у) = иг(1), г = 1, 2, ... п.
(9)
Оно вычислено в тех точках, где было выбрано исходное поле - запись (7). Теперь поля (7) и (9) могут быть сопоставлены между собой. Для минимизации невязок построим такой функционал
Для минимизации функционала используем метод градиентного спуска. Идея этого метода была предложена еще Коши (СоиеИу Огюст Луи 1785-1857). Далее метод совершенствовался и значительное развитие он получил в работах Л.В. Канторовича [1945; 1947; 1948 и др.]. Метод сходится, но точка сходимости всецело зависит от параметров начальной модели - [Булах, 2006]. Для решения задач разведочной геофизики метод получил широкое применение [Страхов, 1964; Булах, 1973; 1976; Старостенко, 1978; Кобрунов, 1995 и др.]. Что касается устойчивости решения, то можно уверенно сказать, что метод устойчив. Достаточно сослаться на работу [Старостенко, Оганесян, 2001].
Все параметры получили начальные численные значения. Начинается итерационный процесс минимизации составленного функционала. Если выполнен к-тый цикл, то (к + 1) получает свои значения
Р
(к + 1) _ (к) л гэ^л
: =Р - Лк и
Коэффициент Лк вычисляется приближенно
Р = Р(№2) = ип(г) - иг(г, №2)]2. (10)
г = 1
Здесь подчеркнуто, что функционал зависит от совокупности параметров. Задача состоит в том, что от их начальных значений W2 = W2(0) необходимо последовательно перейти к параметрам №2*. При этих параметрах невязки сопоставления полей минимизируются.
л Р (Рк) Лк = -
к т Э^
^(эр.
1 = 1 1
Практика показала (см., например [Булах, 1973], что монотонный спуск к решению не гарантируется. Как раз это свойство позволяет улучшить сходимость решения задачи в том случае, ко-
п
Рис. 2. Рельеф дневной поверхности.
гда минимизируется сильно овражный функционал по [Гельфанду и др., 1966].
5. ПРИМЕР
Пусть теперь исследователь располагает исходным полем аномалии силы тяжести - рис. 1. Это поле определено в точках дневной поверхности - рис. 2. На этом рисунке рельеф представлен в изолиниях. Как принято в топографических построениях, аппликата направлена вертикально вверх и возвышенности рельефа имеют большие отметки, чем точки в низинах.
Выберем новую систему координат. Начало этой системы поместим в точке на дневной поверхности в центре участка исследования. Ось аппликат направим вертикально вниз, тогда координатная плоскость XOY будет горизонтальной. Точки аномального
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.