卫星轨道计算：GPS定位过程中的Matlab程序应用

卫星轨道计算在GPS定位过程中卫星轨道位置的定位的前提该程序主要用于对卫星轨道的坐标进行计算采用Matlab编写打开Matlab敲几行代码就能算出头顶这颗GPS卫星此刻在哪这事放在以前得靠天文台蹲望远镜连算好几天对吧毕竟GPS定位说白了就是四个“点光源测距”你不知道卫星自己在哪测出来的距离全是白搭的“瞎数坐标参考线锚点”。卫星轨道计算在GPS定位过程中卫星轨道位置的定位的前提该程序主要用于对卫星轨道的坐标进行计算采用Matlab编写先别着急啃广播星历里那16个飘来飘去的修正项咱们先玩个“纯骨架一点点肉垫”的版本——开普勒轨道地球扁率J2的长期微晃修正这俩凑一块能凑出广播星历90%以上的日常精度够咱们拍个星空定位打卡的虚拟坐标了。第一步先摸清楚这颗卫星的“轨道身份证”轨道根数就是天上的球绕地球转的“入场券动线图说明书”广播星历里每天更新但咱们先挑个真实GPS卫星的典型初值来玩UTC 2024-06-01 00:00:00时刻的PRN01星简化修正后的纯开普勒骨架部分一阶项% 先定义一堆“基础常识常量”别每次查百度了 GM 3.986004418e14; % 地球引力常数GPS官方标准值m³/s² Re 6378137; % WGS84椭球长半轴m J2 1.08263e-3; % 地球扁率二阶带谐系数这个是J2摄动的核心 omega_e 7.2921151467e-5; % 地球自转角速度rad/s % PRN01星的简化轨道身份证骨架J2修正前的 a0 26559700; % 半长轴比广播里多了一点点凑典型 e 0.0012; % 离心率GPS卫星e特别小不然近地点会扎电离层太狠 i0 55.0 * pi/180; % 轨道倾角GPS固定55度±0.1 Omega0 120.0 * pi/180; % 升交点赤经骨架初值 omega0 30.0 * pi/180; % 近地点幅角骨架初值 M0 60.0 * pi/180; % 平近点角初值时刻卫星的平均位置角 t0 0; % 我们把2024-06-01 00:00:00设为t0秒计算方便第二步算J2摄动的“长期小步晃”纯开普勒轨道是固定的椭圆但地球不是完美的球赤道那里鼓了一圈这圈鼓包会拉拽卫星的轨道平面和近地点让Ω和ω每天都转一点点别慌一阶长期项的公式很亲民都是高中三角函数物理量凑比例% 假设我们要算t 3600秒后的位置也就是凌晨1点 t 3600; % 先算平均角速度n0纯开普勒的开普勒第三定律这个熟吧 n0 sqrt(GM / a0^3); % 计算J2导致的Ω和ω的长期变化率dot就是对时间求导的意思 Omega_dot -1.5 * n0 * J2 * (Re / a0)^2 * cos(i0) / (1 - e^2)^2; omega_dot 0.75 * n0 * J2 * (Re / a0)^2 * (5*cos(i0)^2 - 1) / (1 - e^2)^2; % 把骨架初值加上这3600秒的小步晃得到当前时刻的Ω、ω Omega Omega0 Omega_dot * (t - t0); omega omega0 omega_dot * (t - t0); % 轨道倾角i的一阶长期项是0赤道鼓包只会拉平面转不会拉平面歪这个很省心 i i0; % 半长轴a和离心率e的一阶长期项也是0这个椭圆的大小和扁度暂时不变 a a0;第三步算卫星在椭圆轨道平面里的“实时转圈位置”这里要用到开普勒方程——这个方程有点“调皮”是超越方程不能直接用xaxb那种办法解得靠迭代牛顿迭代法最常用收敛贼快3次就够GPS日常精度了% 先算当前时刻的平近点角M就是平均位置角匀速转了t秒 M M0 n0 * (t - t0); % 牛顿迭代解开普勒方程E M e*sin(E) % 初始化迭代GPS的e很小第一次直接用M当E就行 E M; for k 1:3 % 3次足够 f_E E - e*sin(E) - M; % 开普勒方程左边减右边目标让它等于0 f_E_dot 1 - e*cos(E); % 导数 E E - f_E / f_E_dot; % 牛顿迭代公式往0的方向挪一步 end % 看看迭代3次后的E和M差多少肯定小于1e-10 radGPS够用了 fprintf(开普勒方程迭代结果E %.10f radM %.10f rad误差 %.10e rad\n, E, M, abs(E - e*sin(E) - M)); % 算真近点角f这个是卫星和近地点的真实夹角用E转换就行 f 2 * atan2(sqrt(1 e) * sin(E/2), sqrt(1 - e) * cos(E/2));第四步把椭圆平面里的坐标“掰到”地球表面坐标系WGS84地固系椭圆平面的坐标是二维的x轴指向近地点z轴垂直轨道平面咱们要把它转成三维的地固系x轴指向格林尼治子午线z轴指向北极% 第一步算椭圆平面内的二维坐标r是卫星到地心的距离 r a * (1 - e*cos(E)); x_orb r * cos(f); y_orb r * sin(f); z_orb 0; % 第二步三维旋转把轨道平面坐标转成惯性系地心赤道惯性系GCRF不随地球转的 % 旋转顺序是绕z轴转-ω绕x轴转-i绕z轴转-Ω % 旋转矩阵可以分开写也可以合并分开写更清楚 % 绕z轴转-ω的矩阵 Rz_omega [cos(-omega) sin(-omega) 0; -sin(-omega) cos(-omega) 0; 0 0 1]; % 绕x轴转-i的矩阵 Rx_i [1 0 0; 0 cos(-i) sin(-i); 0 -sin(-i) cos(-i)]; % 绕z轴转-Ω的矩阵 Rz_Omega [cos(-Omega) sin(-Omega) 0; -sin(-Omega) cos(-Omega) 0; 0 0 1]; % 合并旋转矩阵注意是从右往左乘先转ω再转i最后转Ω R_GCRF Rz_Omega * Rx_i * Rz_omega; % 算惯性系坐标 pos_GCRF R_GCRF * [x_orb; y_orb; z_orb]; % 第三步把惯性系转成地固系WGS84 ECEF随地球转的 % 从t0到t地球自转的角度是theta_G omega_e * t因为t0设为格林尼治0时刻附近的简化暂时忽略岁差章动春分点这些高阶项够咱们玩了 theta_G omega_e * t; % 绕z轴转theta_G的矩阵 Rz_thetaG [cos(theta_G) sin(theta_G) 0; -sin(theta_G) cos(theta_G) 0; 0 0 1]; % 算地固系坐标 pos_ECEF Rz_thetaG * pos_GCRF; % 输出结果单位换成km更直观 fprintf(PRN01星在UTC 2024-06-01 01:00:00的WGS84地固系坐标\n); fprintf(X %.3f km\n, pos_ECEF(1)/1000); fprintf(Y %.3f km\n, pos_ECEF(2)/1000); fprintf(Z %.3f km\n, pos_ECEF(3)/1000);最后碎碎念两句刚才的代码只加了J2的长期项其实广播星历里还有短期项、周期项、潮汐摄动、日月摄动这些不过日常玩玩虚拟定位或者理解轨道原理这个版本绝对够了。如果想更准直接去IGS官网下高精度的SP3文件然后Matlab里有现成的readsp3函数读就行连自己写迭代都省了。你可以试试把t改成86400秒一天后看看升交点赤经Ω转了多少度——大概-0.04度左右对GPS卫星的轨道平面每天西退大概0.04度这就是为什么你每天同一时间看到的GPS卫星位置稍微偏了一点西边。