Определите положение Сосредоточенного землей зафиксированного землей (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 при помощи Матрицы Направляющего косинуса от 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 при помощи Матрицы Направляющего косинуса от 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: Общая сила тяжести с центробежными эффектами в метрах в секунду придала квадратную форму