CC BY
406
УДК 550.831.01
МЕТОДЫ РЕШЕНИЯ ПРЯМЫХ ЗАДАЧ ГРАВИМЕТРИИ
Ю.В. Пятаков, В.И. Исаев*
Воронежский государственный университет инженерных технологий *Томский политехнический университет E-mail: pyatakovjv@mail.ru
Приведен аналитический обзор методов решения прямых задач гравиметрии. Аппроксимация распределения плотности выполняется набором элементарных тел, имеющих геометрическую форму и заданный закон изменения плотности. Сформулированы требования к оптимальным алгоритмам решения прямых задач и даны рекомендации по их реализации.
Ключевые слова:
Прямые задачи гравиметрии, элементарные тела, закон изменения плотности, оптимальные алгоритмы.
Key words:
Direct gravity problem, the elementary body, the variation of density, optimal algorithms.
Введение
В проблеме развития теории интерпретации геофизических полей важнейшее место принадлежит совершенствованию и развитию методов решения прямых задач.
Решение прямых задач - это определение соответствующих элементов поля по заданным пространственным распределениям физических параметров, характеризующих определенные свойства изучаемой геологической среды (плотность, намагниченность, теплопроводность и др.).
Использование этих методов позволяет учесть влияния геологических тел с известными петро-плотностными и геометрическими параметрами, осуществлять последующий морфологический анализ вычисленных компонент полей, а также решать обратные задачи с использованием методов моделирования.
Задание источников поля должно осуществляться исходя из ряда требований. С одной стороны, оно должно обладать хорошими аппроксима-ционными свойствами, позволяющими наиболее адекватно описывать строение реальных геологических тел и структур. С другой стороны, оно должно обеспечивать высокую скорость вычислений при требуемой точности решения прямой задачи. Требование быстродействия и точности алгоритмов решения прямых задач приобретает первостепенное значение, поскольку их последовательное многократное решение лежит в основе одного из ведущих методов количественной интерпретации - метода подбора. Метод автоматизированного подбора, несмотря на использование регуляри-зованных алгоритмов, может характеризоваться неустойчивостью результатов, что влечет за собой повышенные требования к точности исходных данных и погрешностям вычислений.
В настоящей статье рассматриваются постановки и методы решения прямых задач гравиметрии для типовых аппроксимирующих элементов, а также ряд аспектов, связанных с построением оптимальных алгоритмов их численной реализации.
Постановка прямых задач
Двухмерные задачи гравиметрии. В качестве типового аппроксимирующего элемента в двухмерных задачах удобно использовать трапеции с вертикальными боковыми сторонами и произвольно ориентированными (прямолинейными) верхним и нижним основаниями (рис. 1), а также трапеции, верхнее и нижнее основания которых представлены сегментами парабол (рис. 2).
Рис. 1. Трапеция с вертикальными боковыми сторонами и произвольно расположенными прямолинейными верхним (В) и нижним (Н) основаниями. Условные обозначения и пояснения здесь идалее в тексте
Под прямой задачей гравиметрии для рассматриваемых тел понимается задача определения вертикальной составляющей гравитационного потенциала, выражение для которой имеет вид:
к(„)-2,Др<^>й_ «^*«,
где х и z - координаты точки расчета; f - константа гравитации; £ и С - переменные интегрирования; S - область интегрирования; р(£,£) - распределение плотности.
с(н)(^Мн-Ч2+ь(н)-#+с
н
Рис. 2. Трапеция с параболическими верхним и нижним основаниями
Часто отдельные геологические образования полагают однородными по плотности (т. е. р(£,£)=еош1:), а весь моделируемый объем среды представляется набором кусочно-однородных, различающихся по плотности аппроксимирующих тел.
Вместе с тем, известно, что плотность одновозрастных осадочных пород, вследствие постседи-ментационных процессов, увеличивается с глубиной. Зависимость плотности осадков от глубины наиболее адекватно описывается зависимостями вида [1, 2]:
Р(Я)-В0 -В1ехр(-В,с), (1)
где B0 - минеральная (максимальная) плотность отложений; B1 - характеризует плотность неуплотненных осадков, а B2 - интенсивность воздействия фактора геологического времени.
В случае, когда интервалы изменения глубин отдельных геоблоков невелики, рационально вместо экспоненциального закона (1) использовать его линейную аппроксимацию:
Р(Я)-Ро + К, (2)
где р0 - значение плотности на уровне £=0; k -вертикальный градиент плотности.
Для трапеции на рис. 1 можно полагать, что плотность меняется вдоль вертикального направления по (1) или (2) законам, а для трапеции, представленной на рис. 2 - плотность меняется по закону (2).
Трёхмерные задачи гравиметрии. При решении трехмерных задач в качестве аппроксимирующих тел удобно использовать вертикальные треугольные призмы (рис. 3) с произвольными верхним и нижним основаниями. Плотность в теле может меняться вдоль вертикальной координаты по экспоненциальному (1) или линейному (2) законам.
Выражение для вертикальной составляющей гравитационного потенциала имеет вид:
к (х, у, ^ - / Л/р(С> С-^
в р
где х, у, z - координаты точки расчета; £, ц, £ - переменные интегрирования; D - вертикальная треугольная призма (рис. 3); р(£) - распределение плотности;
Рис. 3. Вертикальная треугольная призма
При построении моделей крупных геологических объектов, когда необходимо учитывать сферичность Земли, в качестве элементов аппроксимации удобно использовать сферические аппроксимирующие тела [3-5]: многогранник (рис. 4) и треугольная призма (рис. 5).
Рис. 4. Сферический многогранник
Определение сферического многогранника. В сферической системе координат £СФ(0;г;р;Я) выделяется тело (рис. 4), имеющее 2xNвершин Я, р. N вершин тела расположены на сфере радиуса г2 и имеют координаты г2, р, Я; /=1,2,..., N а остальные N вершин тела расположены на сфере радиуса г1 (г1<г2) и имеют координаты г1, р, Я; /=1,2,.,N.
Каждая боковая грань тела лежит в плоскости большого круга, проведенного через смежные вершины г2, р, Я; г2, р+1, Я+1 и начало системы координат.
Верхняя и нижняя грани тела образованы соответствующими сегментами сфер радиусов г2 и г1 -сферическими многоугольниками. Каждая сторона (ребро) верхней (нижней) грани представляет собой дугу большого круга, соединяющую смежные вершины г2, р, Я; г2, р+1, Я+1 и гь р, Я; гь р+1, Ям. Боковые ребра тела представлены отрезками прямых, соединяющих вершины с координатами г1, р, Я и г2, р, Я¡. Определенное таким образом тело будем называть сферическим многогранником.
Определение сферической треугольной призмы. В сферической системе координат Есф(0;г;р;Я) выделяется тело (рис. 5), имеющее вершины с координатами (г1в,р1,Я1), (г1н,р1,Я1), (г2в,р2,Я2), (г2н,р2,Я2), (г3”,р3,Я3), (г3н,р3,Я3). Каждая боковая грань тела лежит в плоскости большого круга, проведенного через смежные вершины (гДр,Я); (г,”+1,р+1,Я+1) и начало системы координат.
Рис. 5. Сферическая треугольная призма
Уравнение для верхней грани имеет вид:
г(а)(р, А) - а(а)р + Ь(а)А + с(а).
Соответственно нижняя грань описывается уравнением:
r (i)(p, Я) = a (i)p + b{i)X + c(i).
Плотность выделенных сферических тел полагается постоянной.
Под прямой задачей длярассматриваемых сферических тел понимается определение радиальной составляющей гравитационного потенциала:
vr (д,РоЛ) =
- f рЩ(R - r cos о)P 3r2 sin qdrdtydh. (3)
D
Здесь R, p0, Я0 - координаты точки расчета поля в сферической системе координат; r, р, Я - переменные интегрирования; р - плотность сферического тела D; о- угол при центре сферы между направлениями на точки R, ро, Я0) и г, р, Я, соответственно;
P2 - R2 + r2 - 2Rr cos о, cos о - cos р0 cos р + sin р0 sin р cos^0 - Я).
Алгоритмы решения прямых задач гравиметрии
Создание эффективных методов решения прямых задач - одна из фундаментальных задач теории интерпретации гравитационных полей и ее решению посвящено множество работ.
Идея аппроксимации геологической среды с помощью трапеций с произвольными верхним и нижним основаниями принадлежит Г.Я. Голиздре [6]. В.И. Старостенко и А.Н. Заворотько [7] приводят аналитические выражения для компонент поля от трапеции с постоянной плотностью.
Алгоритмы решения прямых задач для тел с экспоненциальной плотностью можно подразделить на приближенные и точные (аналитические). Приближенные алгоритмы приведены в работах [8—10]. ТВ. Балк с соавторами [8] гравитационный эффект определили путем аппроксимации тела набором полосок, плотность в пределах каждой полоски полагая постоянной и равной значению экспоненты (1) в центре полоски, а гравитационный эффект вычислили как сумму гравитационных эффектов полосок. В [9, 10] тело также разбивалось системой полосок. R. Murthy и B. Rao [9] использовали линейную аппроксимацию экспоненты, позже B. Rao [10] - квадратическую.
В работах [11-13] получены аналитические решения для двухмерных тел (многоугольников [11] и трапеций [12, 13]) с линейной и экспоненциальной плотностью. Ю.В. Пятаковым [12, 13] приведены результаты тестирования и дан анализ полученных выражений в особых точках. Аналитическое решение [12] позволило в явном виде определить выражения для частных производных матрицы Якоби при решения обратной нелинейной задачи. В.И. Старостенко и его учениками [14] выполнена математическая постановка и дано решение обратной нелинейной задачи гравиметрии для сред с экспоненциальной зависимостью плотность-глубина.
В.Н. Страховым [11] получено аналитическое решение прямой задачи гравиметрии для многоугольников с криволинейными сторонами, частным случаем которых является параболическая трапеция (рис. 2).
Аналитические решения прямой трехмерной задачи для многогранников и призм с линейной, а также полиноминальной плотностью приводились в работах [15-17]. Ю.В. Пятаковым и В.И. Исаевым [18] получено и исследовано решение прямой трехмерной задачи для произвольного многогранника с экспоненциально изменяющейся с глубиной плотностью.
Алгоритмы решения прямых задач гравиметрии для сферических тел приведены в работах [3-5]. В.И. Старостенко иА.Г Манукяном [3] получено решение задачи для прямоугольного параллелепипеда, В.Ю. Косыгиным с соавторами [4] - для сферического многогранника, а В.И. Старостенко и О.В. Легостаевой [5] рассмотрен алгоритм решения задачи для сферической призмы с произвольными верхним и нижним основаниями. Интегралы (3) для данных аппроксимирующих тел не выражаются в аналитическом виде. Для их приближенного вычисления использованы квадратурные формулы Гаусса-Лежандра.
Понятие оптимальных алгоритмов
решения прямых задач
Прямые задачи гравиметрии являются классическими и восходят к XIX в. [19]. Для некоторых тел выведены десятки принципиально идентичных, но различающихся по виду вычислительных формул. При этом, одни из формул оказываются чересчур громоздкими и неудобными для программирования на ЭВМ, другие позволяют вычислять поле не во всех внешних по отношению к телу точках, третьи - более универсальны, однако требуют при реализации повышенной точности представления чисел.
В.Н. Страховым [20] были сформулированы основные требования к алгоритмам решения прямых задач: точность, быстродействие и устойчивость. Результаты вычислительных экспериментов [17] позволяют уверенно конкретизировать суть этих требований.
Суть первого требования состоит в том, что вычислительный алгоритм должен обеспечивать изначально заданную точность решения. Поскольку любые вычисления происходят с конечным числом значащих цифр, то существуют ошибки округления. Поэтому, результат, полученный даже по точной аналитической формуле, будет содержать некоторые ошибки, т. е. будет приближенным. В некоторых случаях такие ошибки могут достигать значительных величин и тем самым приводить к неверным результатам. Чтобы избежать этого, необходим анализ ошибок округления и анализ используемой численной схемы.
Требование быстродействия предполагает, что время решения прямой задачи с заданной точностью должно быть минимальным. Несомненно, что максимальную точность и быстродействие удаётся получить в тех случаях, когда для элементов полей оказывается возможным получить явные аналитические выражения в элементарных трансцендентных функциях. Таким образом, проблемы аналитического и точного численного решения прямых задач оказываются взаимосвязанными.
У термина устойчивость для прямой задачи своя специфика. Смысл его раскрывается следующим образом. Существует множество точек пространства, в которых расчеты по аналитическим выражениям без их предварительного анализа невозможны. Вычисляемые величины, входящие в выражения для решения, могут обращаться в бесконечность или претерпевать разрыв, вместе с тем результат должен быть непрерывной функцией. В точках, близких к указанному, так называемому множеству устранимых особых точек, точность вычислений может падать.
Для трехмерных тел в виде однородных многогранников множеством устранимых особых точек является его остов (совокупность ребер и вершин). Аналитическое решение прямых задач гравиметрии для данных тел может быть представлено совокупностью элементарных и трансцендентных функций (логарифмическими и обратными тригонометрическими функциями). Устранимые особенности порождаются логарифмическими функциями, поэтому в работах [19, 21] даны представления их аргументов, позволяющие избежать вычислительной неустойчивости, возникающей при расчетах. В этих же работах рекомендовано использование в качестве обратных тригонометрических функций - функции агС£ и разработана рациональная процедура, обеспечивающая повышение быстродействия вычислительных алгоритмов.
Для сферических трехмерных тел вычислительный алгоритм решения прямой задачи целесообразно строить с использованием квадратурных формул Гаусса-Лежандра, которые обладают наивысшей степенью точности при заданном количестве узлов интегрирования. Для сферического многогранника (рис. 4) интегрирование выражения (3) предлагаем выполнить по двум переменным г, р, что позволяет строить эффективные в плане быстродействия и устойчивости алгоритмы. Сферическая треугольная призма (рис. 5) обладает более универсальными аппроксимационными свойствами (верхняя и нижняя грани этого тела произвольны). Однако аналитическое интегрирование (3) в этом случае можно выполнить только по одной переменной и для численного решения задачи приходится строить кубатурный вычислительный процесс. Для улучшения сходимости вычислительного процесса в точках, лежащих в непосредственной близости к аппроксимирующему телу, предлагаем использовать метод аддитивного выделения особенностей.
Особенности реализации алгоритмов
решений прямых задач
Использование базовых операций FPU. В состав процессоров обычно входит блок операций с плавающей точкой (FPU - Floating-Point Unit), реализующий программно-аппаратную обработку вещественных чисел. Эти устройства в десятки и сотни раз увеличивают скорость выполнения арифметических операций с числами с плавающей точкой. Кроме того, они содержат значения ряда констант и могут выполнять трансцендентные операции.
Оптимальную программную реализацию алгоритмов решения прямых задач предлагаем выполнять с использованием средств языка Ассемблер, поскольку он позволяет осуществить наиболее эффективное применение необходимых базовых команд FPU [22].
Для вычисления обратных тригонометрических функций в алгоритмах предлагаем использовать команду FPATAN [23]. Данная команда имеет два аргумента х и у и вычисляет значение функции arctg (х/у). Результат команды имеет тот же знак, что и операнд х. Время выполнения команды соответствует времени выполнения трех арифметических операций деления.
Для расчета логарифмических функций предлагаем использовать команду FYL2X, которая вычисляет выражение вида: ylog2(x). Пользуясь результатом выполнения этой команды, можно вычислить значение логарифмов по любому основанию.
Помимо операции FYL2X, в блоке FPU имеется команда FYL2XP1, вычисляющая_ выражение ylog2(x+1) при условии: |x|<1-1/V 2. Команда FYL2XP1 обеспечивает высокую точность вычисления логарифмов, когда их аргументы находятся вблизи единицы. В алгоритмах решения прямых задач такая ситуация возникает в тех случаях, когда расчетная точка находится на больших удалениях от аппроксимирующего элемента и возможны потери точности общего решения из-за погрешности вычислений.
Использование технологий параллельного программирования. Расчетные схемы в задачах интерпретации с использованием алгоритмов решений прямых задач обычно организованы по типу вложенных циклических процедур - внешний цикл по точкам расчетной сети системы наблюдений, вложенные циклы - по аппроксимирующим элементам и т. д. При этом одна итерация цикла не зависит от результатов выполнения других итераций. Такие алгоритмические конструкции удобны для реализации в них стандарта OpenMP, позволяющего более полно использовать вычислительные ресурсы многоядерных процессоров [24].
Обсуждение и выводы
Если геологическая структура имеет существенную вытянутость, то можно допустить, что вблизи эпицентра компонента аномального поля, направленная вдоль простирания структуры, отсутствует. К таким структурам относятся ряд рудных и нефтяных объектов. Моделирование гравитационного поля в этом случае целесообразно выполнять с помощью двухмерных аппроксимирующих тел - трапеций.
Если наблюдения проводятся на поисковой площади с линейными размерами порядка десятков и первых сотен км, то направления нормалей к поверхности геоида можно считать практически одинаковыми. Измеренная гравитационная аномалия в этом случае будет представлять собой аномалию вертикальной составляющей гравитационного поля. Для моделирования таких объектов целесообразно использовать вертикальные треугольные призмы.
При решении региональных задач, когда изучаемые структуры имеют протяженность порядка тысячи километров и более, необходимо учитывать сферичность Земли и в качестве аппроксимирующих элементов использовать сферические тела. Сферический многогранник целесообразно использовать для интерпретации данных морских и спутниковых измерений, а сферическую треугольную призму - для данных наземной съемки.
К числу оптимальных алгоритмов, реализующих решение прямых задач для трехмерных тел в виде многогранников, следует отнести отмеченные в настоящей статье алгоритмы В.Н. Страхова. Для сферических тел оптимально использование схемы решения, данное в работах В.И. Старостенко.
Программную реализацию алгоритмов решения прямых задач рекомендуем осуществлять с учетом специфических особенностей архитектуры современных вычислительных средств и технологий программирования. В алгоритмах, реализующих квадратурные вычислительные процессы, рекомендуем использование методов аддитивного выделения особенностей, позволяющих обеспечить устойчивость и быстродействие численных методов.
Авторы статьи занимаются построением и компьютерной реализацией алгоритмов решения задач определения составляющих гравитационного, магнитного, температурного полей, а также распределения скоростей движения и напряжений неоднородных сред. Концептуально принимается, что теория решения трехмерных задач математической физики допускает универсальное использование методов теории потенциалов и многомерных интегральных уравнений для построения алгоритмов эффективного численного решения прямых задач весьма разных по физическому смыслу. Такой концептуальный подход и результаты его реализации будут приведены в следующей статье.
СПИСОК ЛИТЕРАТУРЫ
1. Исаев В.И., Косыгин В.Ю., Лобова ГА., Пятаков Ю.В. Интерпретация данных высокоточной гравиразведки. Вертикальный градиент плотности // Известия Томского политехнического университета. - 2011. - Т 319. - № 1. - С. 83-90.
2. Исаев В.И., Старостенко В.И. Оценка нефтегазоматеринского потенциала осадочных бассейнов Дальневосточного региона поданным гравиметрии и геотермии // Геофизический журнал. - 2004. - Т. 26. - № 2. - С.46-61.
3. Старостенко В.И., Манукян А.Г. Решение прямой задачи гравиметрии на шарообразной Земле // Физика Земли. - 1983. -№ 12. - С. 34-50.
4. Kosygin V.Yu., Pyatakov Yu.V., Krasikov V.N. Direct Gravity Problem for a Spherical Polyhedron // Geology of Pacific Ocean. -1997. - V. 13. - P. 857-868.
5. Старостенко В.И., Легостаева О.В. Прямая задача гравиметрии для неоднородной произвольно усеченной вертикальной прямоугольной призмы // Физика Земли. - 1998. - № 12. -С. 31-44.
6. Голиздра ГЯ. Решение прямой задачи гравиметрии для трехмерных масс // Материалы геофизических исследований на Украине. - Киев: Наукова думка, 1972. - С. 54-66.
7. Старостенко В.И., Заворотько А.Н. Методика и результаты применения регуляризующего алгоритма при решении обратной нелинейной задачи гравиметрии // Геофиз. сб. АН УССР. - 1976. - Вып. 71. - С. 29-40.
8. Балк ТВ., Новоселова М.Р., Балк П.И., Турутанов Е.Х. О точности определения нижней границы геологического объекта по гравиметрическим данным (на примере впадин Байкальской рифтовох зоны) // Физика Земли. - 1988. - № 3. -С. 81-86.
9. Murthy R., Rao B. Gravity anomalies of two-dimensional bodies of irregular cross-section with density contrast varying with depth // Geophysics. - 1979. - V. 44. - № 9. - P. 1525-1530.
10. Rao B. Modelling of sedimentary basins from gravity anomalies with variable density contrast // Geophysics. - 1986. - V. 84. - № 1. -P. 207-212.
11. Страхов В.Н. К теории логарифмического потенциала при переменной плотности возмущающих масс // Физика Земли. -1975. - № 12. - С. 64-81.
12. Пятаков Ю.В. Решение прямой задачи гравиметрии для неоднородных сред // Геофизический журнал. - 1991. - Т. 13. -№5. - С. 79-85.
13. Pyatakov Yu.V., Isaev V.I. Solution of an Inverse Linear 2-D Gravimetric Problem for Sedimentary Cross-Sections // Geophys. J. -1997. - V. 16. - P. 393-406.
14. Старостенко В.И., Пятаков Ю.В., Исаев В.И. Решение обратной задачи гравиметрии для контактов осадочных пород // Геофизический журнал. - 1993. - Т. 15. - № 1. - С. 62-71.
15. Миненко П.А. Моделирование подземных гравитационных аномалий в Кривбассе // Доклады АН УССР Серия Б. -1979. - № 9. - С. 704-709.
16. Кравцов Г.Г. Аналитическое представление внутреннего и внешнего полей притяжения тел с переменной плотностью, ограниченных поверхностями первого порядка // Теория и методика интерпретации гравимагнитных полей. - Киев: Наукова думка, 1981. - С. 105-112.
17. Исаев В.И., Пятаков Ю.В. Решение прямой задачи гравиметрии для трёхмерных блоково-градиентно-слоистых сред // Геофизический журнал. - 1990. - Т. 12. - № 3. - С. 72-79.
18. Пятаков Ю.В., Исаев В.И. Решение прямой задачи гравиметрии для произвольного многогранника с экспоненциально изменяющейся с глубиной плотностью // Геофизический журнал. - 1995. - Т. 17. - № 6. - С. 40-50.
19. Страхов В.Н., Лапина М.И. Прямая и обратная задачи гравиметрии и магнитометрии для произвольных однородных многогранников // Теория и практика интерпретации гравитационных и магнитных полей в СССР. - Киев: Наукова думка, 1983. - С. 3-87.
20. Страхов В.Н. О проблеме решения прямых задач гравиметрии и магнитометрии для материального стержня с полиноми-нальной плотностью // Геофизический журнал. - 1985. -Т. 7. - № 5. - С. 3-9.
21. Страхов В.Н., Лапина М.Н. Прямые задачи гравиметрии и магнитометрии для однородных многогранников // Геофизический журнал. - 1986. - Т. 8. - № 6. - С. 20-31.
22. Магда Ю.С. Использование ассемблера для оптимизации программ на С++. - СПб.: БХВ-Петербург, 2004. - 496 с.
23. Григорьев В.Л. Микропроцессор i486. Книги 2,3,4. - М.: ГРА-НАЛ, 1993. - 382 с.
24. Антонов А.С. Параллельное программирование с использованием технологии OpenMP - М.: Изд-во МГУ, 2009. - 77 с.
Поступила 13.10. 2011 г.
