ДОКЛАДЫ АКАДЕМИИ НАУК, 2013, том 450, № 6, с. 702-707
= ГЕОФИЗИКА
УДК 550.834.3(571.1)
О РЕШЕНИИ ОБРАТНОЙ ЗАДАЧИ ГРАВИМЕТРИИ НА СЕТКАХ БОЛЬШОЙ РАЗМЕРНОСТИ
© 2013 г. Член-корреспондент РАН П. С. Мартышко, И. В. Ладовский, Д. Д. Бызов
Поступило 29.01.2013 г.
БО1: 10.7868/80869565213180187
Для решения задач детального плотностного моделирования необходимо разработать новые компьютерные технологии двух- и трехмерного плотностного моделирования, применимые для обработки числовых массивов большой размерности.
Процесс построения плотностных моделей по аномалиям гравитационного поля сводится к решению прямых и обратных задач гравиметрии. Высокоэффективные алгоритмы "быстрого" решения прямой задачи гравиметрии на сетках большой размерности требуются для успешной реализации функциональных и итерационных схем решения обратных задач. Обратная задача гравиметрии является классическим примером некорректно поставленной: в общем случае ее решение не единственно и неустойчиво зависит от исходных данных. Устойчивое и содержательное решение обратной задачи гравиметрии можно получить, вычисляя отклонения "в малом" от начального распределения плотности модели "нулевого" приближения [1, 2]. Как правило, такие модели строятся по результатам интерпретации профильных сейсмических данных и унаследуют особенности двумерного распределения плотности вдоль сейсмических разрезов. При интерполяции двумерных плотностей в межпрофильное пространство трехмерной модели происходит послойное заполнение сеточных матриц "цифрового куба" [2]. Модель начального приближения, организованная в виде послойных сеточных функций интерполированной плотности, является подходящим объектом для создания устойчивых алгоритмов решения обратной задачи гравиметрии.
Корректное решение трехмерной задачи гравитационного моделирования зависит от способа продолжения масс за боковые грани плотностно-го куба и выбора закона распределения фоновой плотности "нормальной" модели. В качестве та-
Институт геофизики им. Ю.П. Булашевича Уральского отделения Российской Академии наук, Екатеринбург
кой плотности для плоского слоя принимается плотность, зависящая только от глубины. Ее можно получить, усредняя сеточные функции интерполированной плотности по каждому слою "цифрового куба". Такую плотность условно можно назвать гидростатической. Для каждого гипсометрического уровня по плотностной матрице р(х;, у, вычисляются значения средних р0. Полученная совокуп-к
ность средних р0 приравнивается к распределению зависящей только от глубины одномерной плотности р0(х), которая и принимается за плотность относимости "нормальной модели" [2, 3].
ПРЯМАЯ ЗАДАЧА ГРАВИМЕТРИИ ДЛЯ СЕТОЧНОЙ ФУНКЦИИ ПЛОТНОСТИ (ВЫЧИСЛЕНИЕ ЗНАЧЕНИЙ ПОЛЯ ПО ЗАДАННОМУ РАСПРЕДЕЛЕНИЮ ПЛОТНОСТИ)
Введем декартову прямоугольную систему координат, плоскость Оху которой находится на земной поверхности, а ось х направлена вниз. Пусть область Б внутри Земли, заполненная массами с плотностью р(х, у, х), представляет собой прямоугольный параллелепипед, примыкающий к земной поверхности:
Б = ( 0, Ьх )х( 0, Ьу )х( 0, Н);
точка Р(х, у, х) е Б; п, О, С - 0, — точка, в которой вычисляется поле;
Я = ^ - X)2 + (п - у)2 + (С - г)2.
Вертикальная компонента гравитационного поля в точке М вычисляется по формуле (/ — гравитационная постоянная, S — линейный оператор)
g( М) = = йхйуйг. (1)
Б
Воспользуемся сеточной аппроксимацией функции р(х, у, х) плотностного параллелепипеда. На
трехмерной сетке (х, у, гк) построим кубатурные элементы Бку
Б := и, г = 1, 2,..., I, у = 1, 2.....I,
и, к
к = 1, 2.....X,
так, чтобы в пределах каждого из них плотность была постоянной:
р(х, y z ) = pL
_ 1 < X < ,
Уу < У < У г, (2)
£к < г < 1к.
С учетом (2) интеграл (1) заменяется суммой I х /х ^ кубатурных элементов с постоянной плотностью
IX I X X Х Уу
= / X рку | 111^-'^У. (3)
Zk
i,j, к = 1 х-1 y. _ 1
Интегральная сумма (3) в смысле Римана не определена, поскольку существует особенность подынтегральной функции на множестве точек 0 := {R}. Для того, чтобы упростить вычисление вертикальной компоненты гравитационного поля, мы выполним следующие преобразования формулы (3).
Обозначим через Rk расстояние от точки наблюдений M(S, n, Z) до переменной точки интегрирования P(x, y, zk) на глубине zk:
Rk = R\z = zk = V(S - х)2 + (n - y)2 + (Z - Zk)2.
Объединяя в (3) слагаемые с одинаковыми индексами k обратных расстояний, получаем
I х J K
g(S,n,Z) = f X Z , (4)
i, j = 1 к = 0
где Apk. — поэлементная разность плотностей k + 17 а к к + 1 к
и k-го горизонтальных слоев: Api,. = pi,. — pi,.,
k = 1, 2, К- 1, и Ap0. = p1,., ApKj = -pKj ; SkJ -поле вертикального столбца единичной плотности, построенного на ячейках (i, j) с глубины zk до уровня дневной поверхности z0 = 0 [4, 5]:
^ j (S -
х, n - У:
,Z - zk) = i i R - i
dxdy =
-1 yj -1
= {(х - S)ln(y - n + Rk) + (y - n)ln(X - S + Rk) -
- (z - Z) arctg
(х - S ) (y - n ) (z - Z)Rk
(5)
y -1
В первообразной (5) разность логарифмов двух координатных функций следует заменить лога-
рифмом их отношения на двух уровнях глубины. В таком случае элемент поля $ у уже не будет содержать скрытых особенностей на множестве 0 := := [Як], что весьма удобно при вычислениях.
Замечание. Представление решения в форме (4), (5) дает возможность оптимизировать алгоритм вычислений поля для слоя, расположенного на некоторой глубине к ниже уровня дневной поверхности. Для этого достаточно в формуле (5) заменить верхнюю отметку глубины г = 0 на отметку ^ = к и при необходимости сменить индекс к в ку-батурной формуле (4).
Программа вычислений поля для трехмерной сеточной плотностной модели плоского слоя составлена на основе "быстрых" алгоритмов и обладает высокой производительностью. Так, для сеточного куба размером 512 х 512 х 800 узлов время расчета поля на сетке 256 х 256 не превышает 8— 10 мин на стандартном процессоре с тактовой частотой до 3 ГГц.
ОБРАТНАЯ ЗАДАЧА ГРАВИМЕТРИИ ДЛЯ МОДЕЛИ СЛОИСТЫХ СРЕД (ВЫЧИСЛЕНИЕ ПЛОТНОСТИ
ПО ИЗВЕСТНЫМ ЗНАЧЕНИЯМ ПОЛЯ)
Вычисление трехмерной плотности р(х, у, г) неоднородной области Б по заданным на множестве внешних точек значениям поля g(£,, п, С) реализуется на основе решения операторного уравнения первого рода ¿р = g; S — интегральный оператор в правой части формулы (1); g — известная функция. С математической точки зрения такая задача является некорректно поставленной, а ее решение будет сильно зависимым от малых вариаций в исходных данных поля g. Но если на множестве корректности решений обратной задачи выделить класс плотностей, которые меняются только по латерали, то задача о нахождении плотности в горизонтальном слое будет вполне устойчивой [6—8].
Плотность неоднородного параллелепипеда вертикальной мощности Н будем искать в виде произведения зависящей только от глубины функции р0(г) и функции Ф(х, у) [7]:
Р(х, У, г) = Ро(г)Ф(х, у). На разбиении {х;, у, гк} построим сеточный аналог (2) мультипликативной плотности рк,-:
po(z) := {po zk -1 < z < zk}
к = 1'
Ф(х, y) := < Ф, ,■
Xi -1 < х < хi
у. -1 < y < у.
,I х J
Vk !
i,j = 1
х
Поле п, О от слоистого куба на высоте ^ рассчитаем на горизонтальной сетке (г\,у^). На основании (4) имеем
I X J
Щ^О = UWl = X фУ /X APk<j„i,j' (6)
i,j = 1
k = 0
где Apk — разность нормальных плотностей к + 1-и к-го горизонтальных слоев:
AP° = Р° +1- pk, k = 1 2.....К- 1, и
00 Ро = Pk +1 =
Элемент к-го слоя S,j, i;j (поле на сетке (ib j\) от
столбца (i, j) с единичной плотностью) вычисляется по формуле (5). Это означает, что модельный пласт соприкасается с дневной поверхностью. Детальный алгоритм решения линейной обратной задачи гравиметрии будет построен именно для этого частного случая. В силу замечания в конце предыдущего раздела алгоритм автоматически обобщается и на случай погруженного пласта.
Для уменьшения числа индексов оператора (6) прямой задачи введем сплошную нумерацию ячеек (i, j) = X сетки плотностного куба в горизонтальной плоскости:
{i, j }^{X} =: i + (j - 1 ) I,
1 < i < I, 1 < j < J| Л==1X J
В такой же последовательности пронумеруем ячейки (i1, j1) = ц двумерной сеточной функции вычисляемого поля. Размерности сеток {ц} и {X} примем одинаковыми.
Наблюденное гравитационное поле g на дневной поверхности (или его трансформанты на высоте Z) зададим на сетке {ц} расчетного поля U:
g(^il, j Z) = 8». Дискретизация (6) интегрального оператора S прямой задачи на к слоях сетки {X} приводит к матричной форме задачи СЛАУ:
SU Ap0 Фх = . (7)
Здесь Su — трехмерная матрица интегрирования, Apk — приращение плотности по глубине, Ф^ — корректирующая добавка в горизонтальной плоскости.
В уравнении (7) коэффициентами при неизвестных Ф^ являются заранее определенные величины — внутренние суммы по к в правой части формулы (6). По физическому смыслу это значения поля на сетке {ц} вертикальных столбов Dx сетки {X}, набранных из элементов с плотностями pk :
£»,х AP° = S»x = fX APo
(8)
Задача о нахождении плотности в слое сводится к решению СЛАУ вида (7), но меньшей размерности. Двумерная матрица образована сверткой (8) матрицы интегрирования с вектором приращения одномерной плотности и вычисляется только раз. Соответственно, для любого вектора неизвестных Ф^ вычисление поля от трехмерного параллелепипеда сводится к элементарной операции умножения матрицы на вектор:
u» = X .
х = 1
(9)
ИТЕРАЦИОННЫЙ АЛГОРИТМ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ
При разработке устойчивой схемы нами предложен адаптивный итерационный алгоритм решения матричного уравнения (7) линейной обратной задачи, опирающийся на идеи метода локальных поправок [1]. В итерационном методе локальных поправок используется локально-одномерная модель распределения плотности: для каждой итерации приращение сеточного поля от вертикального столба Бх вычисляется только для эпицентральной точки сетки ц = X [1, 6]. Это означает, что из всей суммы (9) для поля учитывается лишь одно слагаемое
(10)
Соответственно, некую эффективную плотность Фх = найденную по эпицентральным значениям поля (10), далее будем называть локально-одномерной плотностью.
Построим итеративный алгоритм после
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.