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

При решении геодезических задач, связанных с определением положения точек на земной поверхности, применяют пространственные геодезические координаты: геодезическую широту B , геодезическую долготу L , геодезическую высоту H . Обработка спутниковых наблюдений при построений глобальных геодезических сетей выполняется в пространственной прямоугольной системе координат X , Y , Z .

Эти системы координат связаны соотношениями [1]:

X = ( N + H ) cos B cos L ;

« Y = ( N + H ) cos B sin L ;

Z = ( N (1 — e2) + H ) sin B;

где N - радиус кривизны первого вертикала; e - первый эксцентриситет земного эллипсоида.

По системе (1) определяются прямоугольные координаты X , Y , Z по заданным геодезическим - B , L , H.

При обратных преобразованиях, т.е. при вычислениях геодезических координат B , L , H по прямоугольным X , Y , Z , наибольшие трудности вызывает определение геодезической широты B при решении трансцендентного уравнения [ 1 ]

tgB = ( Z + e ² N sin B )/ R ; R = V X ² + Y ² .

Для приближенного решения транцендентного уравнения (2), коэффициенты которого являются функциями прямоугольных координат, чаще применяют метод последовательных приближений [1; 2; 3]. Однако при выводе рекуррентных формул многие авторы не уделяют должного внимания, а то и вообще не рассматривают условия сходимости итерационных процессов, скорость приближения их к корню, оценку точности получаемых результатов; не выявляют достоинств и недостатков предлагаемых формул по отношению к известным. Иногда аналитические методы исследования функций подменяются приемом сравнения вычисленных по ним результатов с их более точными величинами для отдельных значений аргументов. Анализ на таком уровне носит поверхностный характер, а выводы и предложения, сделанные на его основе, выглядят малоубедительными, а то и несостоятельными. Для осуществления таких исследований необходимо отделить корень уравнения (2), т.е. указать промежуток как можно меньшей длины, внутри которого он заключен.

С этими же недостатками Госстандартом РФ предложен итеративный алгоритм [4], полученный решением (2) способом введения вспомогательного угла.

Результаты исследований

Для анализа решений уравнения (2) разработана теория и методология отделения корня этого уравнения [5], на основе которой выполнены исследования. Доказано, что уравнение (2) в функциях с приведенной широтой u принимает более простой вид

tgu = (Z1- — e2 + ae2 sin u) Jr .

Эта закономерность позволяет с наименьшими затратами строить алгоритмы и определять решения уравнения (3), а затем полученные результаты с помощью подстановки

tgu = 1- — e2 • tgB

использовать для построения решений уравнения (2).

Из итеративных методов решения уравнения (3) более быструю квадратичную сходимость обеспечивают приближения методом касательных [5]

Z V- — e ² + ae ² sin³ u E tgU =--------2---3---- = ~

.

R — ae cos u     E^

Следовательно, на основе этой математической модели (5) можно построить неитеративное высокоточное решение, представленное только через исходные данные X ,

Y , Z . Для этого величины sin u и cos u , содержащиеся в (5), выразим через X , Y , Z c помощью зависимости [5]

tgu = z 1 - — e ² ■ c₀ Jr ,                                                               (6)

где c о = b J ( b — ae ² ) ;     b = ^R ² + Z ² ( 1 — e ² )                                          (7)

Геометрически величина b – малая полуось воздушного эллипсоида, проходящего через точку Р ( X , Y , Z ) и имеющего тот же эксцентриситет и центр, что и земной эллипсоид.

Для упрощения исследований введем вспомогательные величины:

Д 1 = ^ R ² + ( 1 — e ²) c о Z ² ; Z - = Z /Д - ; R - = R /Д - .                                  (8)

Тогда cos u = 1/ - 1 + tg² u = 1/ ^1 + ( 1 — e ²) c ₀ ² Z²/R ² = R / Д, = R ;

sin u = tgu cos u = c0 V1 — e2 ■ Zx и выражения формулы (5) принимают вид:

E2 = R — ae0 cos3 u = R Ді — ae0 R3 = R (Д! — ae0 R0);

Ex = ZV1 — e0 + ae0 sin3 u = Zx V1 — e0 (Д, + ae0 (1 — e0) c3Z2).

Для дальнейших преобразований E используем зависимость (7), представленную в форме

(1 — e0) c o2 Z 0 = 1 — R2, поэтому

E_x = Z_x 1— — e (Д| + ae c^ — ae c^R_x ) .

С полученными результатами E и E алгоритм (5) принимает вид

Zx V1 — e2  Д, + ae2c0 — ae0c0R2

R '    Д, — ae² R 0

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

Z tgB=R

1+

^c 0 ^{+( 1   c} 0 ) ^R 1

Д 1/ ( ae ²) — R 1 ⁰ _?

Погрешность формулы (9) вычисляем по выражению [5, с. 64]

AB = 3

ae ²

0 ( a + H )

■ sin⁷ B cos³ B .

Исследованиями на экстремум функции (10) устанавливаем, что ее наибольшая величина получается при B = 33°13' и H = 0, составляет AB = 0,0000000" = 2"^10—7, или в линейной мере AB = 0,006 мм. С этими погрешностями AB широта определяется и в точках земной поверхности при |н| < 10 км. С увеличением высоты погрешность AB yбывает и при H ^ ^; AB ^ 0. Для высот спутниковых траекторий H ® 3a получаем AB = 1,9"^10—10.

Для предвычисления точности и расчета вычислительной длины разрядной сетки погрешность (10) выразим через исходные данные с помощью приближенных зависимостей:

a + H = b ₀; sin В = Z/b ₀; cos В = R/b₀ , расчетная формула (10) принимает вид

/

ДВ " =

3(р.р). a-

2      Ч b0 J

( 7

( R

.

.

V b0 J

v b0 J

.

При использовании формулы (11) достаточно вести вычисления с тремя значащими цифрами. Например, по исходным данным:

X = 2773149,800965 м; Y = 1601078,784090 м; Z = 5509233,846979 м получаем b = 636-104; a/b0 = 1,00; Z/b0 = 0,867; Rb0 = 0,504;

Д В = — • 2,76 - 10   - 1 5 • 0,867 7 • 0,504³ = 1,95 - 10 ⁷, или в радианной мере Д В = 10 ¹ рад.

Отсюда следует, чтобы определить широту В по формуле (9) с этой погрешностью, нужно исходные данные X , Y , Z иметь с 13 верными значащими цифрами и вести вычисления, сохраняя один запасной знак. Тогда получим

В = 1,0471975511976 = 60° 00'00", 0000002.

Точное значение В = 60° . Погрешность ДВ = 2" - 10 ^{- 7} подтверждает достоверность как алгоритма (9), так и его оценки (10), (11).

Заключение

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

• 1.    Вычисляем постоянные величины

• 2.    Вычисляем широту В .

к0 = 1 - e2 =(( f -1)/f )2; к = ae2 = a (2 f -1)/ j2, где f - знаменатель сжатия a = (a - b)[a = 1 f.

На эллипсоиде Красовского: a = 6378245 м; f = 298,3.

В геодезической системе координат 2011 г.: a = 637813,5 м; f = 298,2564151.

• 1)    R = X^ 1 + ( YjX ) 2 ;

• 2)    С 1 = к 0 ( Z/R ) ² ;

• 3)    b 0 = R V1 + C 1 ;

• ^{4)    c} 0 = ^b 0/( ^b 0 ^{- k} 1 ) ;

• 5)    Д . = R 7 1 + c . • c 0² ;

• 6)    Z _l = Z / Д₁ ; R _l = R /Д₁ ;

  7) В = arctg

  
  

Z ( 1 + ^c 0 ^{+( 1 -} c 0 ) ^R 1 ² )

R V     Д 1/ k 1 - R ,  J

P.A. Medvedev, M.V. Novorodskaya, S.A. Sharov

Omsk State Agrarian University named after P.A. Stolypin, Omsk

Non-iterative algorithm for calculating geodetic latitude at a spatial rectangular coordinates

When construction of global geodetic networks associated with processing of satellite observations, and the solution of geodetic problems to determine the position of point on the earth surface are applied the spatial rectangular coordinates and geodetic coordinates. In this context arise the problem of transforming these systems of coordinates. This article illustrates the shortcomings of mathematical solutions of the problem of calculating the geodetic latitude B using the method of iterations, having applied the known spatial rectangular coordinates X , Y , Z . Iteration process convergence to an sought quantity does not prove, area and speed of iteration process convergence is not determined, analysis of approaching due to iterations is not carried, comparative study is not carried, advantages and disadvantages of suggested formulas in comparison with existing formulas are not shown. With these drawbacks, the Russian Federation Standard proposes an algorithm for calculating the latitude using the auxiliary angle. On the basis of the previously developed theory of transcendental analysis of solutions a non-iterative algorithm for determination the precise geodetic latitude is derived. To this end, root gap insulation is pre-reduced and a formula for calculating the latitude with high precision by the method of the tangents at one iteration is made. The formula for calculating the error of the proposed algorithm is derived, on the basis of which limits on the accuracy of the results determined by the latitude are set. For the prediction of accuracy and computational calculation of the length of the bit error grid algorithm expressed by the original data. The reliability of the research confirmed in a numerical example.