CC BY
63
УДК 550.831.01
РЕШЕНИЕ ПРЯМЫХ ЗАДАЧ ГРАВИМЕТРИИ ДЛЯ СФЕРИЧЕСКИХ АППРОКСИМИРУЮЩИХ ТЕЛ.
ТЕСТИРОВАНИЕ АЛГОРИТМОВ
В.И. Старостенко, Ю.В. Пятаков*, В.И. Исаев**
Институт геофизики им. С.И. Субботина НАН Украины, г. Киев ‘Воронежский государственный университет инженерных технологий ** Томский политехнический университет E-mail: [email protected]
На системе контрольных примеров выполнено тестирование алгоритмов численного решения прямой задачи гравиметрии для аппроксимирующего тела в виде сферической треугольной призмы с произвольными верхним и нижним основаниями на предмет их устойчивости, точности и быстродействия.
Ключевые слова:
Прямая задача гравиметрии, сферическая треугольная призма, тестирование алгоритмов.
Key words:
Direct gravity problem, spherical triangular prism, testing of algorithms.
Введение
В статье [1] рассмотрены вопросы построения алгоритмов численного решения прямых задач гравиметрии для аппроксимирующих тел, имеющих форму сферического многогранника и сферической треугольной призмы с произвольными верхним и нижним основаниями.
Практическое использование алгоритмов в задачах плотностного моделирования требует предварительного выполнения системы тестовых расчетов для выяснения вопросов, связанных с устойчивостью, точностью и быстродействием данных алгоритмов [2].
В настоящей статье приводятся результаты расчетов и их анализ, выполненные на типовых тестовых примерах, имитирующих работу алгоритмов в условиях их реальной эксплуатации.
Поскольку алгоритм решения прямой задачи гравиметрии для сферической треугольной призмы включает в себя в виде составного блока (см. формулы (11-15) работы [1]) алгоритм решения прямой задачи для сферического многогранника, в статье приводятся результаты тестовых испытаний только для алгоритма решения задачи для сферической призмы.
Тестирование алгоритма решения
прямой задачи гравиметрии
Пример 1. В качестве тестового примера рассматривается сферическая призма, координаты которой по широте и долготе имеют вид: <^=25°, Ф2=25°, ф;=20°, Я1=87,5°, ^=92,5°, Я3=92,5° соответственно (рис. 1). Радиальные координаты вершин верхнего основания призмы (в км): г1в=^+0,2, r2в=R+1,0, rзв=R+1,5. Вершины нижнего основания имеют координаты r1н=R-30,0, r2н=R-32,0 и rзн=R-35,0; R=6371. Верхняя и нижняя грани описываются уравнениями (7) из работы [1] иопира-ются соответственно на вершины верхнего и нижнего оснований призмы. Плотность тела о=1,0 г/см3. Расчетный профиль расположен по
параллели <^=22,5° на расстоянии 1м от верхней сферической грани призмы: га=}"(щ),Л0)+0,001.
В табл. 1 приведены результаты расчета гравитационного эффекта, полученные по квадратурным формулам (11-15) из работы [1], т. е. с использованием метода аддитивного выделения особенностей, и по формулам (10-14) - обычного квадратурного метода, т. е. без выделения особенностей. Точки 1, 2, 3 расположены над телом в непосредственной близости к поверхности его верхнего основания, точки 4, 5, 6 вынесены за проекцию тела на поверхность Земли.
Рис. 1. Проекция сферической треугольной призмы на сферическую поверхность Земли, к тестовым расчетам на примере 1
Количество узлов интерполяции Щ=Щ=П, где п принимает значения 8, 12, 32, 64 и 96. Расчеты осуществлялись на ПК, оснащенном процессором АМБ РИепош II Х4 810 с тактовой частотой ядра 2,8 ГГц.
В точках 4, 5, 6, удаленных от тела, характер сходимости алгоритмов примерно одинаковый. В точках 1, 2, 3, расположенных в непосредственной близости к поверхности тела, как и предпола-
Таблица 1. Результаты расчета радиальной составляющей гравитационного потенциала для примера 1
Расчет по формулам (11-15) из работы [1]
№ точки А Количество узлов интерполяции
8 12 32 64 96
1 90° 9.49762188х102 9,49811659х102 9,49820474х102 9,49820481х102 9,49820502х102
2 91° 1.16546764х103 1,165473 86х103 1,16547386х103 1,165473 86х103 1,16547386х103
3 92° 1,07016303х103 1,07018492х103 1,07018640х103 1,07018640х103 1,07018639х103
4 93° 2,27096899х102 2,27078625х102 2,27077060х102 2,27077060х102 2,27077060х102
5 94° 7,41993108x10' 7,41985079х101 7,41984995х101 7,41984995х101 7,41984995х101
6 95° 3,7047790 5x10' 3,70477365х101 3,70477363х101 3,70477363х101 3,70477363х101
Время счета, с 2,0х10-3 4,3х10-3 О *>< ,9 2, 1,1х10-1 2,5х10-1
Расчет по формулам (10-14) из работы [1]
№ точки А Количество узлов интерполяции
8 12 32 64 96
1 90° 9,4 9723139х102 9,49852010х102 9,49860926х102 9,49857023х102 9,49851431х102
2 91° 1.16542846x103 1,16551425х103 1,16551520х103 1,16551377х103 1,16551132х103
3 92° 1,07013341х103 1,07022523х103 1,07022742х103 1,07022506х103 1,07022123х103
4 93° 2,27083507х102 2,27077373х102 2,27077060х102 2,27077060х102 2,27077060х102
5 94° 7,41986299х101 7,41985003 х10’ 7,41984995х101 7,41984995х101 7,41984995х101
6 95° 3,70477421х10| 3,70477363х101 3,70477363х101 3,70477363х101 3,70477363х101
Время счета, с 1,8х10-3 4,0х10-3 О *>< ,8 2, 1,1х10-1 2,5х10-1
Примечание. Здесь идалее подчеркнуты расхождения срезультатами, выполненными по 96-ти точечной квадратурной вычислительной схеме, реализующей расчет по формулам (1—15) [1]; значения гравитационного поля даны в мгл.
галось в работе [1], сходимость вычислительного алгоритма, реализуемого по формулам (10-14) из [1], очень сильно замедляется.
Проверку точности полученного решения выполним путем аппроксимации исходной призмы системой многогранников. Для этого через начало системы координат в плоскости боковых граней тела проведем прямые, образующие равные углы с прямыми, проходящими вдоль соответствующих боковых ребер призмы. Определим координаты точек пересечения прямых с верхней и нижней гранью исходной призмы и соединим их так, как это показано на рис. 2. Исходное тело разбито таким образом на четыре сферические призмы. Применяя описанную процедуру к каждому телу до тех пор, пока для всех элементов не будет выполнено условие: |ф-ф;-|<£, |А-А7|<£, где /,/=1,2,3, получим разбиение исходной призмы системой элементарных сферических тел.
Рис. 2. Разбиение исходного тела системой элементарных сферических тел
Заменив сферические поверхности элементарных призм плоскими гранями, проходящими через верхние и нижние вершины призм, получим аппроксимацию исходного сферического тела системой многогранников постоянной плотности, гравитационный эффект от которых определим с помощью аналитических формул из работ [3, 4]. Составляющая гравитационного потенциала от системы многогранников рассчитывается в проекции на ось ОТ, связанную с расчетной точкой М.
В табл. 2 приведены результаты расчетов и число аппроксимирующих многогранников для значений £, равных 0,1; 0,01 и 0,001.
Сравнивание результатов расчетов показывают, что использование метода выделения особенностей улучшает сходимость численного алгоритма, не приводя при этом к заметным потерям времени вычислений. Максимальное расхождение результатов расчетов в табл. 1 и 2 - в точке 1, близкой к плоскости боковой грани тела. При 96-иузлах интерполяции алгоритм (11-15) [1] в этой точке обеспечивает гарантировано точность в семь значащих цифр, в остальных точках точность значительно выше.
Пример 2. Координаты тела по долготе и широте, координаты точек расчетного профиля те же, что и в примере 1. Радиальные координаты вершин верхнего основания призмы: г1в=^-30,0;
г2в=^-32,0; г3в=^-35,0. Для нижнего основания г1н=г2н=г3н=420 км. Плотность о=1,0 г/см3.
В табл. 3 приведены результаты расчетов, выполненных по алгоритмам (11-15) и (10-14) [1], соответственно.
Все расчетные точки в данном примере находятся на достаточном удалении от поверхности тела, следствием чего является хорошая сходимость решения, полученного как с помощью алгоритма (11-15), так и с помощью алгоритма (10-14). В табл. 4 приведены
Таблица 2. Результаты расчета гравитационного эффекта в примере 1, полученные путем аппроксимации тела системой многогранников
№ точки Количество многогранников
409б 8б7052 1б77721б
1 90° 9.49796968x102 9,49820572x102 9,49820561x102
2 91° 1,16545458x103 1,16547375x103 1.16547386x103
3 92° 1.07018521x103 1,07018616x103 1,07018632x103
4 93° 2.27088087x102 2,27077105x102 2,27077060x102
5 94° 7,42026297x10' 7.41985174x10' 7,41984991x10'
б 95° 3,70497349x10' 3.70477455x10' 3,70477365x10'
Время счета, с 8,0x10-' 6,6x102 5,1x103
Таблица 3. Результаты расчета радиальной составляющей гравитационного потенциала для примера 2
Расчет по формулам (11-15) из работы [1]
№ точки А:„ Количество узлов интерполяции
8 12 32 64 96
1 90° 2.52184524x103 2.52906160x103 2,52987662x103 2,52987662x103 2,52987662x103
2 91° 2.82394527x103 2,82402742x103 2,82402723x103 2,82402723x103 2,82402723x103
3 92° 2.58088216x103 2.58286480x103 2,58292621 x103 2,58292621x103 2,58292621x103
4 93° 1.94959425x103 1.94760710x103 1,94754571 x103 1,94754571 x103 1,94754571 x103
5 94° 1,40645256x103 1.40644048x103 1,40644052x103 1,40644052x103 1,40644052x103
6 95° 1.04413145x10' 1,04413164x103 1,04413164x103 1,04413164x103 1,04413164x103
Время счета. с 2,0x10-3 4,3x10-3 2,9x10-2 1,1x10"' 2,5x10-'
Расчет по формулам (10-14) из работы [1]
№ точки А:„ Количество узлов интерполяции
8 12 32 64 96
1 90° 2.52191288 x103 2.52906012x103 2,52987662x103 2,52987662x103 2,52987662x103
2 91° 2.82401298x103 2,82402603x103 2,82402723x103 2,82402723x103 2,82402723x103
3 92° 2.58094034x103 2.58286344x103 2,58292621 x103 2,58292621x103 2,58292621x103
4 93° 1.94793039x103 1.94755133x103 1,94754571 x103 1,94754571 x103 1,94754571 x103
5 94° 1,40644093x103 1,40644052x103 1,40644052x103 1,40644052x103 1,40644052x103
6 95° 1,04413160x10' 1,04413164x103 1,04413164x103 1,04413164x103 1,04413164x103
Время счета. с 1,8x10-3 4,0x10-3 2,8x10-2 1,1x10-' 2,5x10-'
Таблица 4. Результаты расчета гравитационного эффекта в примере 2, полученные путем аппроксимации тела системой многогранников
№ точки К Количество многогранников
4096 867052 16777216
1 90° 2.52985442x103 2,52987653x103 2,52987662x103
2 91° 2.82399787x103 2,8240272x103 2,82402723x103
3 92° 2.58290150x103 2,58292612x103 2,58292621 x103
4 93° 1,94753610x103 1,94754567x103 1,94754571 x103
5 94° 1,40643781x103 1,40644051x103 1,40644052x103
6 95° 1.04413098x10' 1,04413163 x103 1,04413164x103
Время счета. с 8,0x10-' 6,6x102 5.1 x103
значения гравитационного поля, определенные с помощью аппроксимации сферического тела системой многогранников. Сопоставление результатов свидетельствует о высокой точности решения.
Пример 3. Известно, что одной из характерных особенностей плотностного строения Земли является градиентное изменение (увеличение) плотности, в зависимости от глубины [5, 6]. Поэтому в данном примере будем полагать плотность тела меняющейся в радиальном направлении. Координаты вершин тела и расчетных точек те же, что и в примере 2. Плотность в теле определим на основа-
нии линейного приближения распределения плотности верхней мантии Земли, полученного по РЕМ-моделям [7] и приведенного в работе [8] в виде соотношения:
a(r) = 7,15855 -3,8599 • r/R (1)
На глубине 35 км значение плотности, определяемое формулой (1), равно 3,32 г/см3, на глубине 420 км значение плотности составляет 3,55 г/см3.
В табл. 5 приведены рассчитанные значения гравитационного эффекта, полученные по алгоритмам (11-15) и (10-14) из работы [1], для тела
Таблица 5. Результаты расчета радиальной составляющей гравитационного потенциала для примера 3
Расчет по формулам (11-15) из работы [1]
№ точки Л Количество узлов интерполяции
8 12 32 64 96
1 90° 8.54924051х103 8.57414809х103 8,57699222х103 8,57699220х103 8,57699220 х103
2 91° 9.56194819х103 9.56223849х103 9,56223779х103 9,56223779х103 9,56223779х103
3 92° 8.74564683х103 8.75254888х103 8,75276503 х103 8,75276503х103 8,75276503х103
4 93° 6.62775163 х103 6.62083436х103 6,62061830х103 6,62061830х103 6,62061830х103
5 94° 4.79844910х103 4,79840637х103 4,79840651 х103 4,79840651х103 4,79840651х103
6 95° 3,57244562х103 3,57244628х103 3,57244628х103 3,57244628х103 3,57244628х103
Время счета, с 2.6 х10-3 5.6х10-3 3,8х10-2 1,5х10-1 3,3х10-1
Расчет по формулам (10-14) из работы [1]
№ точки хч Количество узлов интерполяции
8 12 32 64 96
1 90° 8.54946540х103 8.57414319х103 8,57699222х103 8,57699220х103 8,57699220х103
2 91° 9.56217330х103 9.56223388х103 9,56223779х103 9,56223779х103 9,56223779х103
3 92° 8.74584033х103 8.75254437х103 8,75276503 х103 8,75276503х103 8,75276503х103
4 93° 6.62196251х103 6.62063818х103 6,62061830х103 6,62061830х103 6,62061830х103
5 94° 4.79840799х103 4,79840651 х103 4,79840651 х103 4,79840651х103 4,79840651х103
6 95° 3,57244617х103 3,57244628х103 3,57244628х103 3,57244628х103 3,57244628х103
Время счета, с 2,4х10-3 5,1х10-3 3,7х10-2 1,5х10-1 3,3х10-1
Таблица 6. Сопоставление гравитационных эффектов для тела слинейной и постоянной плотностью
Координата X
90° 91° 92° 93° 94° 95°
V/ 8,5769х103 9,5622х103 8,7527х103 6,6206х103 4,7984х103 3,5724х103
V" 8,6900х103 9.7004х103 8,8722х103 6,6897х103 4,8310х103 3,5865х103
-1,131х102 -1,382х102 -1,195х102 -6,91х10' -3,26х101 -1,41х101
с плотностью, меняющейся в радиальном направлении по линейному закону (1).
В табл. 6 приведены значения гравитационных эффектов V/ и V" соответственно для тела с линейной и постоянной плотностью, а также погрешность ДУ=^УГ'—У"'. При расчете V" плотность полагалась постоянной и равной среднему значению плотности (1) для интервала глубин залегания тела, т. е.
°ср =°(гсР) ~ 3,435 г/см3, гср = Я - (30 + 420)/2 = 6146 км.
По результатам, приведенным в табл. 1, 3, 5, определим среднее время расчета гравитационного эффекта в одной точке для одного тела. В случае постоянной плотности оно составляет соответственно для 96, 64, 32, 12 и 8-и точечной интерполяции, 4,2х10-2, 1,8х10-2, 4,8х 10-3, 7,1х10-4и 3,3х 10-4с. Для тела с плотностью, меняющейся линейно в радиальном направлении, соответствующие значения времени в 1,3 раза больше. Как видно из результатов расчетов, минимальное количество интерполяционных узлов, необходимых для достижения заданной точности, зависит от расположения точки по отношению к телу (по мере удаления от тела число узлов уменьшается).
Пример 4. На рис. 3 приведен сферический треугольник с размерами 1°х1°. Радиальные координаты точек тела те же, что и в примере 1, координаты по широте и долготе ф;=ф2=21°; ф3=20°; Я;=89,5°; ^=90,5°; Я3=90,0°. Верхняя и нижняя грани призмы
описываются уравнениями (7) из работы [1] и опираются соответственно на вершины верхнего и нижнего оснований. Плотность тела о=1,0 г/см3. Расчетные точки расположены на разных удалениях от тела, так как это показано на рис. 3. Все точки находятся на расстоянии одного метра от поверхности верхней грани призмы, т. е. г0;=гв(ф0, Ао;)+0,001.
Рис. 3. Проекция призмы и расчетных точек на сферическую поверхность Земли в примере 4
Таблица 7. Результаты расчета радиальной составляющей гравитационного потенциала для примера 4
№ точки Количество узлов интерполяции
8 12 32 64 96
1 5,67089975х102 5,67087331х102 5,67087373x102 5,67087376x102 5,67087363x102
2 1,80880724х102 1,80882432х102 1,80882403 x1°2 1,80882403 x1°2 1,80882403 x1°2
3 6,34505817x10' 6,34505795x10' 6,34505795x10' 6,34505795x10' 6,34505795x10'
Время счета, с 1,0х10-3 2,1 х10-3 1,4 x1с-2 0,5x10-' 1,2x10-'
№ точки Количество узлов интерполяции
4 6 8 12 32
4 2.81702998x10' 2,81702578x10’ 2,81702550x10' 2,81702550x10' 2,81702550x10'
5 2.80579960x10° 2,80580152х100 2,80580152x100 2,80580152x100 2,80580152x100
6 2,78660749x10° 2,78660929x100 2,78660930x100 2,78660930x100 2,78660930x100
7 2,58938055х100 2,58938152x100 2,58938152x100 2,58938152x100 2,58938152x100
8 2,1298899х100 2,12989033x100 2,12989033x100 2,12989033x100 2,12989033x100
9 5,60378376х100 5,60376221x'0° 5,60376209x100 5,60376209x100 5,60376209x100
10 5,43614749х100 5,43612937x100 5,43612928x100 5,43612928x100 5,43612928x100
11 4,84659990х100 4,84658458x100 4,84658452x100 4,84658452x100 4,84658452x100
12 4,09333253х100 4,09332388x100 4,09332386x100 4,09332386x100 4,09332386x100
Время счета, с 3,1х10-3 6,4x10^ 4,3 x10-2 1,6x10-' 3,6x10-'
Результаты расчетов данного примера, приведенные в табл. 7, позволяют определить степень зависимости минимального необходимого числа узлов интерполяции, обеспечивающих заданную точность, от положения точки по отношению к телу.
Из табл. 7 видно, что только для точек 1, 2 и 3, расположенных в непосредственной близости к телу, потребовалось более 8-и узлов интерполяции. Для остальных точек шесть и восемь узлов обеспечивают достаточно высокую точность (относительные погрешности составляют величины порядка 10-6и10-8, соответственно).
При моделировании объектов, представленных большим количеством аппроксимирующих элементов (несколько сотен и более) основной объем вычислений будет приходиться на расчет поля в точках, удаленных от соответствующих тел. Поэтому при оценке времени вычислений представляется возможным ориентироваться на среднее время вычислений, выполненных по 8-и точечной схеме интерполяции. Например, расчет гравитационного поля модели слоя из 1000 аппроксимирующих тел в
3000 точках должен потребовать не более 17 мин.
процессорного времени для указанного выше ПК.
Выводы
1. Разработанные алгоритмы решения прямых задач гравиметрии для принятых сферических тел обеспечивают высокую точность результата не зависимо от места расположения точки расчета по отношению к аномалиеобразующему объекту.
2. Использование метода аддитивного выделения особенностей при расчетах, осуществляемых в непосредственной близости точки расчета к поверхности тела, существенно повышает точность результатов по сравнению с обычным квадратурным методом.
3. Для точек расчета, удаленных от аппроксимирующего элемента, количество интерполяционных узлов, обеспечивающих заданную точность, закономерно уменьшается. Достижима высокая скорость расчетов при обработке данных реальных сложнопостроенных структур (с использованием тысяч аппроксимирующих элементов).
СПИСОК ЛИТЕРАТУРЫ
1. Старостенко В.И., Пятаков Ю.В. Решение прямой задачи гравиметрии для сферических аппроксимирующих тел. Алгоритмы // Известия Томского политехнического университета. -2013. - Т 322. - №1. - С. 28-34.
2. Пятаков Ю.В., Исаев В.И. Методы решения прямых задач гравиметрии // Известия Томского политехнического университета. - 2012. - Т. 320. - № 1. - С. 105-110.
3. Страхов В.Н., Лапина М.И. Прямая и обратная задачи гравиметрии и магнитометрии для произвольных однородных многогранников // Теория и практика интерпретации гравитационных и магнитных полей в СССР. - Киев: Наукова думка, 1983. - С. 3-87.
4. Страхов В.Н., Лапина М.Н. Прямые задачи гравиметрии и магнитометрии для однородных многогранников // Геофизический журнал. - 1986. - Т. 8. - № 6. - С. 20-31.
5. Старостенко В.И., Легостаева О.В. Прямая задача гравиметрии для неоднородной произвольно усеченной вертикальной прямоугольной призмы // Физика Земли. - 1998. - № 12. -С. 31-44.
6. Исаев В.И., Косыгин В.Ю., Лобова Г.А., Пятаков Ю.В. Интерпретация данных высокоточной гравиразведки. Вертикальный градиент плотности // Известия Томского политехнического университета. - 2012. - Т. 319. - № 1. - С. 83-90.
7. Жарков В.Н., Трубицин В.П. Физика планетных недр. - М: Наука, 1980. - 448 с.
8. Хачай Ю.В. Исследования конвекции сжимаемой гравитирующей жидкости в плоском слое применительно к условиям верхней мантии Земли // Методы интерпретации и моделирования геофизических полей. - Свердловск: УрО АН СССР, 1988. - С. 82-89.
Поступила 03.09.2012 г.
