Определение положения Earth-Centered Earth-Fixed (ECEF)

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

Вычислите геоцентрические радиусы в метрах в массиве широт от +/- 90 градусов с помощью geocradius.

lat = -90:90;
r = geocradius( lat );
rlat = convang( lat, 'deg', 'rad');
z = r.*sin(rlat);

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

x = r.*cos(rlat);
y = zeros(size(x));

Вычислить зональную гармоническую гравитацию для Земли

Использовать gravityzonal вычислить массив зональной гармонической гравитации в координатах ECEF для массива позиций ECEF в метрах в секунду в квадрате.

[gx_zonal, gy_zonal, gz_zonal] = gravityzonal([x' y' z']);

Расчет силы тяжести WGS84

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

lat_gd = geoc2geod(lat, r);
long_gd = zeros(size(lat));
[gd_wgs84, gn_wgs84] = gravitywgs84( zeros(size(lat)), lat_gd, long_gd, 'Exact', [false true false 0] );

Определение силы тяжести для сферической земли

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

GM  = 3.986004415e14;
gd_sphere = -GM./(r.*r);

Графики сравнения для различных гравитационных моделей

Для сравнения гравитационных моделей их выходы должны находиться в одной системе координат. Можно преобразовать зональную гравитацию из координат ECEF в координаты NED, используя Матрицу направления косинуса (Direction Cosine Matrix) из dcmecef2ned.

dcm_ef = dcmecef2ned( lat_gd, long_gd );
gxyz_zonal = reshape([gx_zonal gy_zonal gz_zonal]', [3 1 181]);
for m = length(lat):-1:1
    gned_zonal(m,:) = (dcm_ef(:,:,m)*gxyz_zonal(:,:,m))';
end

figure(1)
plot(lat_gd, gned_zonal(:,3), lat, -gd_sphere, '--', lat_gd, gd_wgs84, '-.')
legend('Zonal Harmonic', 'Spherical', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Down gravity (m/s^2)')
grid

Рис. 1: Сила тяжести по оси вниз в метрах в секунду в квадрате

figure(2)
plot( lat_gd, gned_zonal(:,1), lat_gd, gn_wgs84, 'r-.' )
hold on
plot([lat(1) lat(end)], [0 0], 'g--' )
legend('Zonal Harmonic', 'WGS84', 'Spherical','Location','Best')
xlabel('Geodetic Latitude (degrees)')
ylabel('North gravity (m/s^2)')
grid
hold off

Рис. 2: Гравитация по оси Север в метрах в секунду в квадрате

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

gtotal_wgs84 = gravitywgs84( zeros(size(lat)), lat_gd, zeros(size(lat)), 'Exact', [false true false 0] );
gtotal_zonal = sqrt(sum([gx_zonal.^2 gy_zonal.^2 gz_zonal.^2],2));

figure(3)
plot( lat, gtotal_zonal, lat_gd, gtotal_wgs84, '-.' )
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Total gravity (m/s^2)')
grid

Рис. 3: Общая гравитация в метрах в секунду в квадрате

Сравнение гравитационных моделей с центробежными эффектами

Вы видели гравитационные сравнения невращающейся Земли. Изучите центробежные воздействия вращения Земли на модели гравитации.

Вычислить гравитационные центробежные эффекты для Земли

Использовать gravitycentrifugal вычислить массив центробежных эффектов в координатах ECEF для массива положений ECEF в метрах в секунду в квадрате.

[gx_cent, gy_cent, gz_cent] = gravitycentrifugal([x' y' z']);

Добавление центробежных эффектов к зональной гармонической гравитации.

gx_cent_zonal = gx_zonal + gx_cent;
gy_cent_zonal = gy_zonal + gy_cent;
gz_cent_zonal = gz_zonal + gz_cent;

Расчет WGS84 гравитации с центробежными эффектами

Использовать gravitywgs84 для вычисления WGS84 гравитации по нисходящей и северной осям на поверхности Земли массив геодезических широт в градусах и 0 градусов долготы с использованием точного метода с атмосферой, центробежными эффектами и без предварительной обработки.

[gd_cent_wgs84, gn_cent_wgs84] = gravitywgs84( zeros(size(lat)), lat_gd, long_gd, 'Exact', [false false false 0] );

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

gtotal_cent_wgs84 = gravitywgs84( zeros(size(lat)), lat_gd, zeros(size(lat)), 'Exact', [false false false 0] );
gtotal_cent_zonal = sqrt(sum([gx_cent_zonal.^2 gy_cent_zonal.^2 gz_cent_zonal.^2],2));

Графики сравнения для различных гравитационных моделей с центробежными эффектами

Для сравнения гравитационных моделей их выходы должны находиться в одной системе координат. Можно преобразовать зональную гравитацию из координат ECEF в координаты NED, используя Матрицу направления косинуса (Direction Cosine Matrix) из dcmecef2ned. На рисунке 5 видно некоторое различие между зональной гармонической гравитацией с центробежными эффектами и WGS84 гравитацией с центробежными эффектами. Большая часть разницы обусловлена различиями между зональной гармонической гравитацией и WGS84 расчётами гравитации.

gxyz_cent_zonal = reshape([gx_cent_zonal gy_cent_zonal gz_cent_zonal]', [3 1 181]);
for m = length(lat):-1:1
    gned_cent_zonal(m,:) = (dcm_ef(:,:,m)*gxyz_cent_zonal(:,:,m))';
end

figure(4)
plot(lat_gd, gned_cent_zonal(:,3), lat_gd, gd_cent_wgs84, '-.')
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Down gravity (m/s^2)')
grid

Рис. 4: Гравитация с центробежными эффектами по оси вниз в метрах в секунду в квадрате

figure(5)
plot( lat_gd, gned_cent_zonal(:,1), lat_gd, gn_cent_wgs84, '--', lat_gd, (gned_zonal(:,1)-gn_wgs84'), '-.' )
axis([-100 100 -0.0002 0.0002])
legend('Zonal Harmonic', 'WGS84', 'Error Between Models w/o Centrifugal Effects', 'Location','Best')
xlabel('Geodetic Latitude (degrees)')
ylabel('North gravity (m/s^2)')
grid

Рис. 5: Гравитация по оси Север в метрах в секунду в квадрате

figure(6)
plot( lat, gtotal_cent_zonal, lat_gd, gtotal_cent_wgs84, '-.' )
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Total gravity (m/s^2)')
grid

Рис. 6: Общая гравитация с центробежными эффектами в метрах в секунду в квадрате