Этот пример показывает, как исследовать зональное гармоническое, сферическое и 1 984 Мировых Геодезических Системы (WGS84) модели силы тяжести для широт от +/-90 градусов в поверхности Земли.
Поскольку система координат 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']);
Используйте 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;
Используйте 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: Общая сила тяжести с центробежными эффектами в метрах в секунду придала квадратную форму