Приближенные методы глобального гармонического сферического анализа потенциальных полей Текст научной статьи по специальности «Математика»

УДК 550.831

И. В. Бойков, М. В. Кравченко

ПРИБЛИЖЕННЫЕ МЕТОДЫ ГЛОБАЛЬНОГО ГАРМОНИЧЕСКОГО СФЕРИЧЕСКОГО АНАЛИЗА ПОТЕНЦИАЛЬНЫХ ПОЛЕЙ

Аннотация. Предложен алгоритм определения коэффициентов Фурье в разложении по сферическим функциям потенциальных полей в предположении, что потенциальные поля известны своими значениями на неравномерной сетке узлов, заданной в некоторой части поверхности Земли. Метод основан на использовании полиномов Бернштейна для продолжения потенциальных полей на всю поверхность Земли и последующем применении кубатурных формул. Ключевые слова: коэффициенты Фурье, сферические функции, потенциальные поля, полиномы Бернштейна, кубатурные формулы.

Abstract. In the paper the algorithm of definition of Fourier coefficients in decomposition of potential fields on spherical functions. The potential fields are known by its values on a non-uniform grid of units, given in some part of the Earth surface.

The method is based on the Bernshtein multinomials using for continuation of potential fields on all surface of the Earth and subsequent application of the cubature formulas.

Keywords: Fourier coefficients, spherical functions, potential fields, Bernshtein polynomials, cubature formulas.

Введение

Задачи глобального гармонического сферического анализа и синтеза восходят к классическим работам К. Гаусса [1, 2] и Ф. Нейманна [3, 4] по теории земного магнетизма.

Несмотря на то, что и К. Гаусс, и Ф. Нейманн проводили вычисления вручную, предложенные ими алгоритмы представляют интерес и в настоящее время.

И К. Гаусс, и Ф. Нейманн при построении глобального сферического гармонического синтеза использовали двухступенчатый метод.

Алгоритм К. Гаусса заключался в том, что на первом шаге применялось преобразование Фурье по переменной ф (здесь используется сферическая система координат (р,0,ф), в которой 0 - долгота, ф - широта). На втором шаге - метод наименьших квадратов.

Алгоритм Ф. Нейманна отличался от алгоритма К. Гаусса на втором шаге. Вместо метода наименьших квадратов он применял квадратурные формулы наивысшей алгебраической степени точности.

Метод Нейманна получил в дальнейшем широкое распространение в различных разделах геофизики. Он используется в топографии [5, 6], в физике атмосферы [7], в геодезии [8-10].

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

Приведем краткий обзор работ по двуступенчатому методу в задачах глобального сферического гармонического синтеза (ГСГС) и глобального сферического гармонического анализа (ГСГА).

Пусть функция f (0, ф) на сферической поверхности S разлагается в ряд по сферическим функциям

то n

f (0, ф) = 22 Р^ (cos 0)(Cnm cos тФ + Snm sin тф) =

n=0m=0

TO TO

= 22 K(n, m)РЩ1 (cos 0)(Cnm cos тф + Snm sin тф) =

n=0m=0

TO TO

= 2 2 Pm (cos 0)(Cnm cos 1ф + Snm sin 1ф) =

m=0n=m

TO

= 2 ^m (0)cos 1ф + Bm (0)sin 1ф, (1)

m=0

TOTO

где 4m (0)= 2 Pnm (cos 0)(Cnm ), Bm (0) = 2 P^ (cos 0)(Snm ), K(n, m) = 1 при

n=m n=m

n > m, K(n, m) = 0 при n < m, Pnm (cos 0) - нормированные присоединенные функции Лежандра,

рт /<2л + 1) (n-m)! (1 -12)m/2 dn+m(t2 - 1)n

n \ 2 (n + m)! 2nn! dtn+m ’

где t = cos0.

Глобальный сферический гармонический анализ заключается в нахождении коэффициентов Cnm и Snm из разложения (1) функции f (0, ф) по сферическим функциям.

Используя ортогональность тригонометрических полиномов и присоединенных полиномов Лежандра, из разложения (1) имеем

\Cnm\ =--------1--------ff f (0, ф)Pm (cos 0) Icos 1ф I ds, (2)

ISnm J (1 + §m0)^JSJ 1 sin 1ф\ ’ ^

где ds = sin0d 0d ф, S = [0,:rc;0,27i], 8m0 = -

5

1,т = 0 [0,т Ф 0

Применяя к формуле (2) кубатурные формулы, можно вычислить коэффициенты Спт и 5пт.

Двухступенчатый ГСГА состоит в том, что последовательно проводятся вычисления по переменным ф и 0. В двухступенчатом ГСГС вычисления проводятся в обратном порядке - сначала проводятся вычисления по 0, а затем - по ф.

Остановимся на этих алгоритмах подробнее.

Алгоритм двухступенчатого ГСГА заключается в том [11], что вначале вычисляются функции

ґ ~\ 2п г

\Лт (0) ] 1 1 f ч[С08 тф

1 j f (0, ф) {С08гаф1 ¿Ф, (3)

п J [ sin тф J

[Bm (0) J (1 + §т0)--0

а затем вычисляются коэффициенты

(^il^1±^i(Bm(e)}im(cose)si"erfe- (4)

Двухступенчатый ГСГС заключается [11] в том, что вначале проводится суммирование по параметру n

М-р'НЙ- (5)

а затем вычисляются значения функции

то

/(в, ф) = ^ (Am (0)cos тф + Bm (0)sin тф). (6)

т=0

Алгоритмы сферического анализа и синтеза развивались в работах [9,

12, 13].

Непосредственное применение формул (3), (4) возможно только при небольших значениях n и m. Поэтому необходима дискретизация этих формул.

Основная проблема при этом заключается в том, что после дискретизации присоединенные полиномы Лежандра оказываются неортогональными.

Это накладывает дополнительные трудности на решение задачи ГГСА даже при условии, что дана равномерная сетка узлов на всей сферической поверхности.

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

Для вычисления отрезков ряда по шаровым функциям

1 Nnpnm (cos 0), . .

U(Х) = ~ ----П---(anm C0S ^ + bnm sin

г n=0m=0 г

где v>0 - целое число; P^1 (cos0) - присоединенные функции Лежандра; x = (х1,х2,хз), (г,ф,0) - сферические координаты, в [14] построены четырехчленные рекурсивные соотношения, связывающие присоединенные функции Лежандра Pn, нормированные следующим образом:

Pnm (cos 0) = -^ Pnm (cos 0).

(n + m)!

Алгоритм «столбцового» типа позволяет проводить устойчивый синтез потенциальных полей при N > 180 с высокой точностью [14].

В работе [15] описана методика синтеза потенциальных полей на персональных компьютерах в среде MathLab™, позволяющая осуществлять синтез до N = 75.

Как отмечалось выше, в ГСГА потенциальных полей широко используется метод наименьших квадратов. Однако этот метод, как показано в [16], требует очень большого числа наблюдений.

Метод наименьших квадратов эффективно применяется в случае, если задана равномерная по широте и долготе сетка наблюдений. Однако в геодезической практике такая сетка отсутствует, в частности, из-за проблемы полярных областей. В работе [13] отмечается, что в полярных областях отсутствует регулярная сетка спутниковых наблюдений и ее приходится дополнять наземными и авиационными съемками.

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

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

В данной работе предложен метод решения этой проблемы.

1. Вспомогательные предложения

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

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

Рассмотрим сегмент [0,1], на котором имеется N узлов

0 = хо < x <... < Xn =1. Пусть Ak = [xk, xk+l], k =0,1,..., N -1, hk = xk+1 - xk,

*

k = 0,1,.., N -1. Будем считать выполненным условие h /h* < m, где

*

h =maxhk, h*=minhk, k = 0,1,..,N -1, m - целое число.

Построим следующую модификацию полиномов Бернштейна:

N

Pn (х) = 24/(Xk) xk (1 - х)N-k. (7)

k=0

Повторяя рассуждения, приведенные в [18], можно показать, что для непрерывной функции /(х), х е[0,1], lim ||/(х) - Bn (х)|| ^ 0 при выполне-

*

нии условия h /h* < m, m = const.

Точно также, повторяя доказательство теоремы Т. Поповичиу [18, с. 245246], можно показать, что для модифицированных полиномов Бернштейна справедлива оценка

\/ (х) - ВN (х) \< 3т / ;^=

где ю (/, 8) - модуль непрерывности функции / (х).

Известна [18, с. 254] теорема Л. В. Канторовича, утверждающая, что если / (х) есть целая функция, то ее полином Бернштейна

Вм (х) = ¿4/( N ) хк (1 - х)м-к к=0 ^ м'

сходится к ней на всей числовой оси.

Повторяя рассуждения, приведенные в [18, с. 254-256], можно показать, что если и(х) - сужение на числовую ось гармонической функции, то последовательность полиномов Бернштейна сходится к этой функции.

2. Численный двухступенчатый алгоритм ГГСА, основанный на экстраполяции полей и применении кубатурных формул

Введем сферическую систему координат (г, 0, ф) с центром в центре сферы S радиуса Я. Пусть известны значения потенциального поля и (Я, 0, ф) в области О на поверхности сферы, определяемой неравенствами

0< «1 <0<«2<л, 0<Р1 <ф<Р2 < 2л. Положим Р1 = 0,Р2 = 2л, т.е.

0 <ф< 2л. Ограничение, связанное с предположением, что 0 <ф< 2л, не влияет на общность рассуждений.

Предлагаемый алгоритм ГСГА состоит из двух этапов.

Первый этап

На первом этапе, располагая значениями потенциального поля на поверхности , продолжаем его на поверхность сферы S радиуса Я .

Пусть значения функции и(Я, 0, ф) заданы на прямоугольной сетке

узлов {0к,ф/} , к = 1,2,...М1, I = 1,2,...,N1, причем эта сетка неравномерная

ни по переменной 0, ни по переменной ф . По узлам {0к,ф/} , к = 1,2,...,М1,

/ = 1,2,...,N1, построим полином Бернштейна:

Вад^, ф) = 2 2 сМ1 сМ1и ( я, 0к, ф/)

к=1 / =1

0 - «1 I I 1 0 - «1

М1 М | 0 « |к | 0 « \М1—к

,сы, см1и (

( ф 0 I1( ф 0 I м1-/

а 2 -«1

а 2 -а1

х

х

°2 - °1 ) I °2 - °1

Полином Бернштейна Вм1 м1 (0, ф) аппроксимирует функцию и(Я, 0, ф) в области О.

Теперь введем в области й равномерную сетку узлов:

К +(а2 _а1 )-----5 Р1 + (Р2 _Р0- Г 5 k, 1 = П1 •

[ П1 П1 )

По значениям полинома Бернштейна Вм (0, ф) на равномерной сетке узлов строится новый полином Бернштейна:

Вщт(0, ф) - 2 2 СП СП Дц,^ а1 + (а2 - а1)—, рі + (р2 - Рі)—

1 1 1 1 1 1 ^ Пі Пі

—-01-0

X

X

0- аі а2 - аі

Ґ

і-

0- аі а2 - аі

\Пі-— ґ а Ґ

Ф-Рі Р2 -рі

і-

Ф- Рі Р2 -Рі

\Пі-1

Замечание. Необходимость в построении этого полинома обусловлена тем, что теорема Канторовича о сходимости последовательности полиномов Бернштейна к целым функциям доказана для полиномов Бернштейна, построенных на равномерных сетках узлов.

*

Полином ВП1 п (0, ф) аппроксимирует функцию и(Я, 0, ф) на сфере 51.

Замечание. В случае, если нужная точность восстановления функции

*

и(Я, 0, ф) полиномом Вп1 П1 (0, ф) не достигается во всей области S одновременно, то, как отмечается в [17], можно построить последовательность полиномов Бернштейна, определенных на последовательности областей, сходящейся к S.

*

Построением полинома Вп п (0, ф) заканчивается первый этап алгоритма.

Второй этап

*

Располагая значениями полинома Вп п (0, ф) на поверхности S, вычислим коэффициенты Фурье по сферическим функциям. Для этого воспользуемся формулой

К

|5П

і

(і + 8 Ф)РТ(С080)Р^ IД.

(і + 8,ио)я" 1,1 {вт тф \

(8)

Перейдем в (8) к сферической системе координат. Воспользовавшись формулой интегрирования по сфере [і9], имеем

С

К

і

п 2 л

И* - т Геов тф!

В?2і ,Пі (0, Ф)РП (сов 0) 1 . Г БШ м0^ф.

^■~т0/'“00 і {в1П тФ і

(9)

Так как значения функции Вп^ п1(0, ф) легко вычисляются при любых значениях (0, ф), то для вычисления интеграла из правой части формулы (9) можно использовать различные кубатурные формулы.

і06

По переменной ф естественно воспользоваться квадратурной формулой по равноотстоящим узлам, которая является формулой наивысшей тригонометрической степени точности. Для вычисления интеграла по переменной 0 можно воспользоваться квадратурной формулой Гаусса п -го порядка.

В результате получаем формулу

ІС

Ч— 17

1

2л

^пш \ (1 М ^_і і_і

тіі^гіі

X РПЩ I 008 X

(N) к

2пш1

008-

81П-

М

2пш1

М

8ІП ХкМ) + Яш (и).

(10)

.(N) х(NД

где ', - коэффициенты и узлы квадратурной формулы Гаусса;

РЫМ (и) - погрешность квадратурной формулы (10).

При реализации кубатурной формулы (10) возможны два подхода. Во-первых, можно взять достаточно большое значение М, а значение N взять равным п . При этом при каждом значении п приходится брать новые наборы весов и узлов квадратурной формулы Гаусса. Во-вторых, можно взять достаточно большое значение N и М .

Наряду с формулой (10) вычисление коэффициентов Фурье сферических функций можно проводить по кубатурным формулам с равноотстоящими узлами:

Сп

К

1 2л к к* (пк 2яТ| х

(1+§шо)~МмЬ\к п1,п11N.М х

X РЩ\ 008

пк

N

2пш1

008

81П-

М

2пш1

М

пк

81П N + ККМ (и) •

(11)

3. Модельный пример

Введем декартову систему координат ОХХХ. Пусть в точке с координатами (30,30,30) находится точечное тело с массой ц = 100. создающее потенциальное поле и (х. у. г).

Введем сферическую систему координат (р. 0. ф) с полюсом в начале декартовой системы координат. Пусть потенциальное поле и (х. у. г) известно на поверхности {р = Я = 10.10°<0< 80°. 0 <ф< 2л} и равно

и (х. у. г) =

100

д/(X - 30)2 + (у - 30)2 + (г - 30)2

(12)

Требуется определить коэффициенты Фурье сферических функций (2). Для этого будем использовать численный алгоритм, приведенный во втором разделе. Полученные результаты сравниваются с точными значениями коэффициентов Фурье разложения поля по сферическим функциям:

2(2 — 8гао)Ц R * к m/ А*ч

Cnm =----------—-—cosтф Pn (cos 0 ),

2n +1

/ *\ n+1 (r )

Snm = ^ 8_0)^ R„, , sinтф*рпт(cos0 ),

2n +1

¡ *\n+1 (r )

(13)

* * * s

где (г , 0 , ф ) - координаты источника потенциального поля. Результаты вычислений приведены на рис. 1.

n

Рис. 1. График погрешности коэффициентов Фурье На рис. 1 p(n) = max(| Cnm — Cnm |), где Cnm - коэффициент, найден-

m

ный по формуле (10); Cnm - точное значение коэффициента, найденное по формуле (13). Сплошной линией изображен график функции p(n) при использовании равномерной сетки узлов для построения полиномов Бернштейна. Штрихпунктирной линией показан график функции p(n) при использовании неравномерной сетки узлов для построения полиномов Бернштейна *

при соотношении h /h* < 2 .

С помощью найденных коэффициентов Cnm и Snm восстановим функцию U (R, 0, ф) по формуле для внешней краевой задачи [20]:

N' n / r \ n+1

UN' (r, 0, ф) = 2 2[ —I Pnm (cos 0) (nm cos mф + Snm sin mф), r > R .(14)

n=0 m=0 ^ Г '

В табл. 1 представлена максимальная погрешность восстановления функции U(R, 0, ф), заданной формулой (12), по формуле (14) при использо-

вании равномерной и неравномерной сетки узлов для построения полиномов Бернштейна.

Таблица 1

N ' Формула Равномерная сетка Неравномерная сетка

Максимум относительной погрешности, % Максимум абсолютной погрешности Максимум относительной погрешности, % Максимум абсолютной погрешности

10 Гаусса (10) 2,49746E-001 4,28310E-003 3,58967E-001 7,48856E-003

прямоугольников (11) 1,12281E+000 2,34652E-002 1,08009E+000 2,06412E-002

20 Гаусса (10) 2,21755E-001 4,60688E-003 3,62166E-001 7,59521E-003

прямоугольников (11) 3,99273E+000 8,43455E-002 3,91479E+000 7,63034E-002

30 Гаусса (10) 2,31276E-001 4,85023E-003 3,47737E-001 7,29261-003

прямоугольников (11) 8,92589E+000 1,90004E-001 8,88127E+000 1,79664E-001

Выводы

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

Список литературы

1. Гаусс, К. Ф. Избранные труды по земному магнетизму / К. Ф. Гаусс. - Л. : Академия Наук СССР, 1952. - 344 с.

2. Гаусс, К. Ф. Избранные геодезические сочинения / К. Ф. Гаусс. - М. : Геодезиздат, 1957. - 144 с.

3. Neumann, F. Uber eine nene Eigenschaft der Laplaceschen у(") und ihre Anwendung zur analytischen Darstellung derjenigen Phänomene / F. Neumann // Schumachers Astron. Nachr. - 1838. - V. 15. - P. 313-325.

4. Neumann, F. Vorlesungen uber die Theorie des Potentials und der Kugelfunctionen / F. Neumann. - Leipzig : Teubner, 1887. - P. 135-154.

5. Prey, A. Darstellung der Hohen- und Tiefenverhaltnisse der Erde durch eine Entwicklung nach Kugelfunktionen / A. Prey // Math. Phys. Kl. Neue Folge. - 1922. -V. 11. - № 1. - P. 134-167.

6. Hofsommer, D. J. On the Expansion of a Function in a Series of Spherical Harmonics / D. J. Hofsommer. - Amsterdam : Computation Department of the Mathematical Centre, 1957. - 344 p.

7. Ellsaesser, H. W. Expansion of Hemispheric Meteorological Data in Antisymmetric Surface Spherical Harmonic (Laplace) Series / H. W. Ellsaesser // J. Appl. Meteorology. - 1966. - № 5. - P. 263-276.

8. Payne, M. H. Truncation Effects in Geopotential Modelling / M. H. Payne. -Maryland: Analytical Mechanics Associates, 1971. - 367 p.

9. Colombo, O. L. Numerical Methods for Harmonic Analysis on the Sphere / O. L. Colombo. - Ohio : The Ohio State University. Department of Geodetic Science and Surveying, 1981. - 310 p.

10. Пеллинен, Л. П. Высшая геодезия (теоретическая геодезия) / Л. П. Пеллинен. -М. : Недра, 1978. - 264 с.

11. Sneeuw, N. Global spherical harmonic analysis by least squares and numerical quadrature methods in historical perspective / N. Sneeuw // Geophysical Journal International. - 1994. - V. 118. - № 3. - P. 709-716.

12. Rizos, C. An Efficient Computer Technique for the Evaluation of Geopotential from Spherical Harmonic Models / C. Rizos // Aust. J. Geod. Photo. Surv. - 1979. - V. 31. -P. 161-169.

13. Sanso, F. Fast spherical collocation theory and examples / F. Sanso, C. C. Tscherning // Journal of Geodesy. - 2003. - V. 77. - P. 101-112.

14. Страхов, В. Н. Методы синтеза рядов по шаровым функциям, представляющих элементы потенциальных геофизических полей / В. Н. Страхов, А. Б. Ефимов, М. М. Хохрякова // Известия АН СССР. Физика Земли. - 1988. - № 5. -С. 41-57.

15. Bethencourt, A. Using personal computers in spherical harmonic synthesis of high degree earth geopotential models / A. Bethencourt, J. Wang, C. Rizos, A. H. W. Kearsley // Dynamic Planet. - 2005. - P. 125-130.

16. Мориц, Г. Современная физическая геодезия / Г. Мориц. - М. : Недра, 1983. -392 с.

17. Бойков, И. В. О приближенном методе восстановления потенциальных полей / И. В. Бойков, М. В. Кравченко, В. И. Крючко // Известия РАН. Физика Земли. -2010. - Т. 46. - № 4. - С. 67-77.

18. Натансон, И. П. Конструктивная теория функций / И. П. Натансон. - М. : Государственное издательство технико-теоретической литературы, 1949. - 688 с.

19. Крылов, В. И. Приближенное вычисление интегралов / В. И. Крылов. - М. : Наука, 1967. - 500 с.

20. Тихонов, А. Н. Уравнения математической физики / А. Н. Тихонов, А. А. Самарский. - М. : Наука, 2004. - 798 с.

Бойков Илья Владимирович

доктор физико-математических наук, профессор, заведующий кафедрой высшей и прикладной математики, Пензенский государственный университет

E-mail: [email protected]

Кравченко Марина Витальевна

аспирант, Пензенский государственный университет

E-mail: [email protected]

Boykov Ilya Vladimirovich Doctor of physico-mathematical sciences, professor, head of sub-department of higher and applied mathematics,

Penza State University

Kravchenko Marina Vitalyevna Postgraduate student,

Penza State University

УДК 550.831 Бойков, И. В.

Приближенные методы глобального гармонического сферического анализа потенциальных полей / И. В. Бойков, М. В. Кравченко // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 4 (16). - С. 101-110.