Решение прямых задач гравиметрии для сферических аппроксимирующих тел. Алгоритмы Текст научной статьи по специальности «Физика»

Геофизика

УДК 550.831.01

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

АЛГОРИТМЫ

В.И. Старостенко, Ю.В. Пятаков*

Институт геофизики им. С.И. Субботина НАН Украины, г. Киев *Воронежский государственный университет инженерных технологий E-mail: pyatakovjv@mail.ru

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

Ключевые слова:

Прямая задача гравиметрии, сферический многогранник, сферическая треугольная призма.

Key words:

Direct gravity problem, spherical polyhedron, spherical triangular prism.

Введение

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

При построении моделей крупных геологических объектов необходимо учитывать сферичность Земли, а расчеты проводить в сферической системе координат [1-3]. В качестве элементов аппроксимации удобно использовать сферические аппроксимирующие тела [1-4]. В работах [1-3] для аппроксимации плотностного строения среды используется сферический прямоугольный параллелепипед с плотностью, меняющейся в радиальном направлении по линейному закону, в работе [4] плотностная параметризация среды осуществляется с использованием сферической прямоугольной призмы с произвольно расположенными верхним и нижним основаниями. Плотность аппроксимирующего элемента в работе [4] меняется вдоль параллелей и меридианов пропорционально длинам дуг на верхнем и нижнем основаниях, линейно или экспоненциально вдоль любого радиуса и принимает заданные значения в вершинах призмы.

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

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

Прямая задача гравиметрии

для сферического многогранника

Постановка задачи. В соответствие с работой [5], дадим определение рассматриваемому аппроксимирующему элементу.

В геоцентрической системе координат Есф(0;г;фД) (определение которой приводится в работах [1-3]) выделим тело Б (рис. 1, а), имеющее 2xN вершин.

N вершин тела расположены на сфере радиуса г2 и имеют координаты тъ ф, X /=1,2,..,N остальные N вершин расположены на сфере радиуса г1 (^!<г2) и имеют координаты гь ф, X.

Каждая боковая грань тела лежит в плоскости большого круга, проведенного через смежные вершины г2, ф, X; г2, Фи, Х+1 и начало системы координат.

Верхняя и нижняя грани тела образованы сегментами сфер радиусов г2 и гх- т. н. сферическими многоугольниками. Каждая сторона (ребро) верхней (нижней) грани представляет собой дугу боль-

шого круга, соединяющую смежные вершины г2, ф, X; r2, ф;+1, Xi+1 и г1, ф, X; г1, ф;+1, Я;+1. Боковые ребра тела представлены отрезками прямых, соединяющих вершины с координатами г1, ф, X и г2, ф, X. Определенное таким образом тело будем называть сферическим многогранником.

Под прямой задачей для сферического многогранника D будем понимать определение радиальной составляющей гравитационного потенциала:

у, (я,ф„А) =

= f fflp(r)(R - г cos а)P Зг2 sin (pdrdydX. (1)

D

Здесь R, ф0, X, - координаты точки расчета поля в сферической системе координат; D - сферический многогранник; г, ф, X - переменные интегрирования; а - угол при центре сферы между направлениями на точки R, ф0, X, и г, ф, X соответственно;

P2 = Я2 + г2 - 2Rг cos а, cos а = cosф0 cosф + sin ф0 sin фcos(X0 -X).

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

р(г) = р + к(,2 - г), (2)

где р - постоянная плотность; к - коэффициент, характеризующий изменение плотности по направлению г.

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

теграл в правой части (1) можно представить в виде суммы интегралов:

f fflpir)( R -r cos а)P 3r2 sin<pdrd<pd)i =

D

N

=E ‘i• f JJJpir)(R - r cos а)P 3r2 sinфdrdфdЯ, (3)

i=l D,

где

‘ i = sign[iaix0 + b Уо + Ciz0)iaixi+2 + b іУ i+2 + c izi+2)];

a = y,zi+l-y,+lzi; bt = -xizi+l + x,+lz,; c, = x y *+l- w,; x, =sin ф cos Я; yt =sin ф sin Я;

zi = cosФi; i = l,2,...N+2; ф^і =Фі; ФN+2 =Ф2;

Я+l = Л; -^N+2 = Я2.

Расчет интегралов в правой части (З) удобно проводить в системе координат OX'YZ', полученной ортогональным вращением системы координат OXYZ вокруг центра O таким образом, чтобы точка M оказалось на оси OZ', как это показано на рис. 2, а.

Тогда выражения для слагаемых в правой части (З) можно записать в виде:

f fflpir )і R -r cos а) P 3r2 sin фdrdфdЯ =

D,

= f fflpir)(R -r cos ф')P 3r2 sin фdrdфdЯ', (4)

Di

где Dt - элементарный сферический многогранник (рис. 1, б, 2, а); ф', Я' - переменные интегрирования, определенные в системе координат ЕСФ(0;г;ф';Я), связанной с декартовой системой координат OX'YZ'.

Подставляя в правую часть (4) выражение (2) для плотности р(г) и выполняя последовательное интегрирование по переменным ф' и г, получим выражение следующего вида:

/ Щр(г)(К - г соБф') Р V 2Бт ф'drdф'd X' = ц

= /(<г + к ■ г2)К|^(XVX -/ ■ к ■ К21Г2(ХУX’, (5) 0 0

где

F1(X’) =

= t

(2 - 3t)Pk К2Pt +t • rlPk = 2 + 1

3 3 k = 1

+2t (1 -12)Arth A(t)

t = cosф'(1 ')

t = 1 ’

F2(X' ) =

=t

2

.J

3t Pk (rk -t) + 3rkPk + 5t • Pk 2 4 4

3(1 -12) pk (r'-1)

-rk3 Pk

k = 2

t = СОБф ' (1 ') t=1;

к = 1

+3(1 - / 2)2ЛгШ Л( /)/4

гк = гк/К; Рк =л11+К2 - 2/1 ■ к ;

Л(?) = (г2 - г^/( Р2+Р1).

Воспользовавшись формулами сферической тригонометрии [6], определим выражение для углов Ц (рис. 2, а):

Ц = 2arctg

[cos(at - bt) - cos ct ] [cos c t - cos(at + bt)]

где

cosat = cosф0cosф +sinф0sinф cos( 10-1); cosbt = cosф0cosф+1 +sinp0sinф+1cos(10 -1+1); cos ct = cos ф cos ф+1 +sinф sinф.+1^( 1 -1+1). Аналогично определим выражение для cos^'(1):

cos ф' (1') = (1 - рф2) /^(1 + Д2)(1 + Р2),

где

р, = sin[(A -1 ) / 2] • tg(bt / 2) / sin[( A + 1')/2].

Р2 = cos [(Л -1) / 2] • tg(bt / 2) / cos[( A +1) /2], tg (bt /2) = 7(1 - cos bt )(1 + cos bt). sin[( A ±1' )/2] =

= sin(A/2)^cos(1/2) ±cos(A/2) ^sin( 1/2). sin( A /2) =*J[cos(bt - ct) - cos a | / 2 / sin b / sin ct. cos(A /2) = ^/[cos at -cos(bt + ct)]/ 2 / sin b / sin c. cos[( A ±1')/2] =

= cos(A/2) ^cos(1/2) + sin( A/2) ^sin( 1/2).

Л - угол при г-й вершине сферического треугольника (рис. 2, б).

Интегралы в правой части (5) не выражаются в конечном виде, поэтому для их вычисления будем использовать квадратурную формулу Гаусса-Лежандра, обладающую наивысшей алгебраической степенью точности при заданном количестве узлов интегрирования [1-4, 7, 8]:

Lt т n

J Fk 1 )d1 * 2E AjFk (1j). k = 1.2; t = 1.2..... N. (6)

0 2 j=1

где 1j=Lx/2; n - порядок квадратуры; xи Aj - узлы и коэффициенты квадратурной формулы.

Для обеспечения расчета интегралов по формуле (6) с заданной точностью используем процедуру, описанную в работе [3].

Рис. 2. К расчету радиальной составляющей гравитационного потенциала для элементарного сферического многогранника D,. а) определение угла Ц; б) определение угла ф'(Х')

При к=0 решение (1-6) дает, как частный случай, решение задачи для сферического многогранника постоянной плотности, рассматриваемой в работе [9].

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

Постановка задачи. Как и в случае сферического многогранника, в соответствии с работой [5], дадим определение рассматриваемому аппроксимирующему элементу.

В геоцентрической системе координат Есф(0;г;ф;Х) выделим тело (рис. 3), имеющее вершины с координатами (г1в,ф1,Х1), (г1н,ф1,Х1),

(г2в,Ф2,Х), (г2н,Ф2,Х), (Гзв,Фз,Хз), (гзн,Фз,Хз). Каждая боковая грань тела лежит в плоскости большого круга, проведенного через смежные вершины г1в, ф, X; гв+1, ф+1, Х+1('=1,2,з; Г4в=Г1в; Г4н=Г1н; ф4в=ф1в; Х4в=Х1в) и начало системы координат.

Верхнее и нижнее основания призмы (поверхности) опираются на верхние и нижние вершины тела и описываются уравнениями:

гв (ф, X) = а вф + ЬВХ + св;

г н(ф, X) = а нф+ Ь НХ+ с н. (7)

Поверхности (основания) Г‘(ф,Х), г"(ф,Х) будем называть сферическими, а рассматриваемое тело -треугольной сферической призмой с произвольно расположенными верхним и нижним основаниями.

Будем полагать, что плотность призмы меняется в радиальном направлении по линейному закону:

р(г) = р- к ■ г. (8)

Под прямой задачей гравиметрии для сферической треугольной призмы Б будем понимать опре-

деление радиальной составляющей гравитационного потенциала:

Vr (^ф,А) =

= f JJJp(r)(R -rcosa)P 3r2sinqdrdqd 1. (9)

D

Здесь R, ф0, Xi, r, ф, 1 - те же, что и в выражении (1); p(r) - плотность, описываемая соотношением (8).

Аналогом данного аппроксимирующего элемента в декартовой системе координат является вертикальная треугольная призма с плотностью, меняющейся с глубиной по линейному закону -типовой аппроксимирующий элемент, рассмотренный в работе [10].

Решение задачи. Интеграл в выражении для гравитационного эффекта сферической призмы допускает аналитическое интегрирование только по одной переменной, поэтому для численного решения прямой задачи (9) будем использовать квадратурный вычислительный метод, описанный в работах [1-4]. Так же как и при решении прямой задачи для сферического многогранника, выполним поворот системы координат так, чтобы расчетная точка M оказалось на оси OZ новой системы координат. Подставим в правую часть (9) выражение для плотности (8) и, выполнив интегрирование по переменной г, получим следующее выражение для Vr(R^o,1o):

V r (Л.фоД,,) = f JJ F (ф ',1')sinф' dф' d1' .

(10)

где

= p[R • F32 -cosф' • F33] -k[R • F33 -cosф' • F34]. (11)

r (ф>1) rm

Fm = Fm (ф 1) = J —dr;

n n Vr > / J pn

т=0,1,2,з,...; и=2/+1; г=0,1,2,_ (12)

Выражения для интегралов (12) записываются с помощью рекуррентных соотношений, приведенных в работе [з]:

т--0,5 тт-1р_2-1ьрт-1 -

m - 2t

m-1

m 1 rm-1 p-2-

m - 2t

F2 = —b/2F2'~‘+F^2 -r 1 2t+1 ' 2t+1 ~ 1 2-1 1

-/P-

2+1

m - 2t

2 -1 / (2 t +1;

F =

-1 2t+1

2(2r+b)

2t -1 8(t -1)

/ AP

-21+1

/ (2t A* 0.

(2t + 2)-1 P

A = 0;

F0 =-ln(P — b/2-r); a =R2; b =-Rcosф; A = 4a - b2.

В зависимости от расположения расчетной точки М относительно тела, показанного на рис. 4, расчетная формула (10) может принимать одну из следующих форм:

3 1т ф'т (Х )

а) Гг(К,ФоА) = /^ \ | F(ф\X')sinф’dф^Х

т=1 о о

1+1 1т ф'т (Х' )

б) Гг(К’оХо) = /| ^(ф'Х>тФ^Ф^X';

т=1 о о

1} ф} ^ ')

в) Гг(Яф^о) = /1 | F(ф',X')sinф’dф'dX';

о о

1] ф}(X')

г) Гг(R,ф0,X0) = /1 | F(ф\X')sinф'dф'dX'.

0 ф} (X')

1] ф} (X')

д) Гг(Я,ф0Д0) = /1 | F(ф',X})sinф^ф’dX +

0 ф} (X')

^ фк (X ')

+/$ | F^ф’,x ')sin ф ' dф ' dX X .

о ф}(X')

Здесь Ц, фj'(X) - те же, что и в формуле (5) для сферического многогранника.

Значения f(ф,X), гн(фД) в выражениях (12) определяются соотношениями (7), представленными в исходной системе координат ZСФ(0;г;ф;X). Для определения этих значений в системе координат SСФ(0;г;ф';X) приведем формулы расчета для ф(ф',X), X(ф'X), полученные с помощью формул [7]:

ф(ф' ,X' ) = arccos d, X(ф',X' ) = аг^(5/ с), (1з)

где

а б

d = cosф' cosф0 + sin ф sin ф0cos 1

5 = sinф' / sinф ' • sin 1 .

c = (cos ф' - cos ф0 cos ф^/sin ф0 / sin ф' .

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

1 ф2

вида J J F (ф ' .1')sin ф' dq'd1'. который будем осу-

0 фi

ществлять с помощью квадратурных формул Гаусса-Лежандра:

1 ф2

J J F(ф'.1 ')sin ф 'dф 'd 1' *

0 ф

* 1 /4 • E (ф (1,') - ф (1,0)£ AA F(fp'j .1 1) sin ф'. (14)

t=1 j=1

где

1,' = 1( xt +1)/2;

ф' = (ф2 (1/ ) -ф1 (1 ))х, / 2 + (ф2(1 ' ) + ф(1 ' )) / 2;

X, х,- - узлы; Аг, Aj - весовые коэффициенты квадратурных формул Гаусса-Лежандра.

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

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

в

Рис. 4. Положение точки относительно тела: а) над телом; б) над телом в плоскости боковой грани; в) над телом на продолжении бокового ребра; г, д) за проекцией тела на поверхность Земли (в случае г точка принадлежит плоскости боковой грани тела). І, к - номера вершин

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

Подинтегральная функция в (10) содержит слагаемые вида г") 8тф'/Г, где Р - расстояние от точки М до точки поверхности основания призмы с координатами

г(ф', X' ),ф' , X' ;

Р = л]г2(',X') + К2 -2соъ'-К ■ г',X'); г(ф',X') = гв(<р,X); ф=ф(ф',X'); Х=Х(<р',Х').

При ф'^0 величина Р^0. Поскольку 8Шф' и Р имеют одинаковый порядок малости, особенность г&пф'/Р является устранимой, однако, учитывая повышенные требования метода Гаусса-Лежандра к точности вычисления подинтегральных функций, а также к их гладкости, желательно исключить из расчетов точки, в которых соответствующие функции могут иметь особенности. Для этого удобно воспользоваться известным в вычислительной математике методом аддитивного выделения особенностей [10].

Пусть расчетная точка М (Я,'^!) находится вблизи поверхности г"") (рис. 5) или непосредственно на самой поверхности. Вычислим значения Г1=1‘ (фоЛ), г=1* " и построим сферический многогранник Б 'координаты вершин которого ф, (/=1,2,3) совпадают с соответствующими

координатами призмы Б, а радиальные координаты равны г1 и г2 соответственно.

3

\ Л/(7?;ф0;/Ч1)

Рис. 5. 1) сферическая призма D; 2) сферический многогранник D'; 3) поверхность верхнего основания призмы

Преобразуем формулу (10) следующим образом:

Кг ( к'X =

= ^И [^(''XX) - ^' (ф'X)^Шф 'ёф'dXX +

Б

+У'г (К,ф0,X0), (15)

где ^'^ф^!) - гравитационный эффект многогранника Б, вычисленный по формулам (1-6), а

Г(ф',Х) вычисляются аналогично Р'(ф',Х) по формулам (11), (12), в которых вместо г" (ф0,Х0) положено r2, а вместо гн(ф,Х)-г1. Расчет интеграла в (14) выполняем с помощью квадратурной формулы (14). Соответствующие слагаемые в подинтеграль-ном выражении в (15) будут иметь вид si^'(ry_P-r 2/Р ), где P - расстояние от расчетной точки М до точки с координатами г(ф',Х), ф', Х', лежащей на поверхности верхней грани сферической призмы Б, а P'- расстояние от Мдо соответствующей точки, принадлежащей верхней грани сферического многогранника Б' (рис. 5). При фф——0, r—r, а P—P, таким образом, соответствующие слагаемые в подинтегральной функции в (15) при ф— 0 стремятся к нулю, что должно обеспечить хорошую сходимость метода квадратур.

Теперь рассмотрим алгоритм решения прямой задачи для случая, когда плотность тела может изменяться как по переменной r, так и по ф и Х, принимая заданные значения в вершинах призмы:

р(ф,Х,г) = р(ф,Х) - к(ф,Х) • г, (16)

где

р(ф,Х) = р0 + р ф + р2 ‘Х, к(ф, Х) = к0 + kj ф + к2 •Х, (17)

р0, р, р2, к0, кь к2 вычисляются с помощью соответствующих интерполяционных формул [8].

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

vr (Д,ф„А) =

= f Шр(ф, Х, r)( R - г cos a)P 3r2 sin фdrdфdХ. (18)

D

где р(ф0,Х,г) определяется соотношением (16).

Численный алгоритм решения задачи (18) для плотности (16, 17) аналогичен алгоритму решения задачи (9) для плотности (8). Выражение (11) для функции Дф',Х) здесь будет иметь вид:

F(ф',Х') = р(ф,Х)[R • F2 - ^ф' • F33] -

-к(ф,Х)[R • F33 - cosф'• F34],

где зависимости ф=ф(ф',Х), Х=Х(ф',Х) определяются формулами (13).

В формуле (15) Vr'(Я,ф0,Х0) - гравитационный эффект сферического многогранника с плотностью, меняющейся вдоль радиального направления по линейному закону р(г)=р(ф0,Х0)-к(ф0,Х0)г.

Заключение

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

2. Разработан алгоритм решения прямой задачи гравиметрии для аппроксимирующего тела в виде сферического многогранника с плотно-

стью, меняющейся в радиальном направлении по линейному закону.

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

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

4. Алгоритмы решения задач реализуются с использованием численных методов: квадратурного метода Гаусса-Лежандра и (в прямой задаче для сферической призмы) метода аддитивного выделения особенностей.

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

СПИСОК ЛИТЕРАТУРЫ

1. Старостенко В.И., Манукян А.Г Решение прямой задачи гравиметрии на шарообразной Земле // Известия АН СССР. Фи-

зика Земли. - 1983. - № 12. - С. 34-50.

2. Старостенко В.И., Манукян А.Г. Задачи гравиметрии для луче-

вых ускорений и изучение маскона Моря Нектара на Луне // Известия АН СССР. Физика Земли. - 1988. - № 7. - С. 60-72.

3. Старостенко В.И., Манукян А.Г., Заворотько А.Н. Методика

решения прямых задач гравиметрии и магнитометрии на шарообразных планетах - Киев: Наукова думка, 1986. - 112 с.

4. Старостенко В.И., Легостаева О.В. Прямая задача гравиметрии для неоднородной произвольно усеченной вертикальной прямоугольной призмы // Известия АН СССР. Физика Земли. -1998. - № 12. - С. 31-44.

5. Пятаков Ю.В., Исаев В.И. Методы решения прямых задач гравиметрии // Известия Томского политехнического университета. - 2012. - Т. 320. - № 1. - С. 105-110.

6. Кранц П. Сферическая тригонометрия - М: ЛКИ, 2007. - 98 с.

7. Страхов В.Н. Метод приближенного решения прямой трехмерной задачи гравиметрии // Известия АН СССР. Физика Земли. - 1983. - № 9. - С. 52-62.

8. Бахвалов Н.С., Жидков Н.П., Кобельков ГМ. Численные методы. - М: БИНОМ, 2006. - 636 с.

9. Kosugin V.Yu., Pyatakov Yu.V., Krasikov V.N. Direct Gravity Problem for a Spherical Polyhedron // Geol. of Pac. Ocean. - 1997. -V. 13. - P. 857-868.

10. Пятаков Ю.В., Исаев В.И., Косыгин В.Ю. Методы теории потенциала при решении прямых задач гравиметрии и геодинамики трехмерных неоднородных сред // Известия Томского политехнического университета. - 2012. - Т. 321. - № 1. -С. 76-83.

Поступила 03.09.2012 г.