《天际回响：中继卫星轨道确定》（含美国中继卫星系统(TDRSS)系统仿真代码）

美国的跟踪和数据中继卫星系统（Tracking and Data Relay Satellite,TDRSS）为所有主要的美国空间观测站、研究卫星以及载人航天飞机提供跟踪服务。尽管该系统独具特色，但它为卫星-卫星跟踪技术提供了一个典型的示例。在这里，我们将讨论跨多个转发器的信号路径建模以及多航天器轨道调整的方法。

1 数学模型

如原理图所示，中继卫星的测距测量是通过发送测距信号到一个地球同步中继卫星开始的。从这里，它被转发到目标用户航天器，再次重发后通过中继卫星回到地面站。设定在地面站接收信号的时间为 ，用户卫星、中继卫星和地面站的惯性位置向量分别为 、 和 ，则可以得到如下隐含方程：

通过这些关系从初始值为零开始逐次确定光程时间。得到的值则提供了模型化的四路测距：

忽略地面站和两颗卫星的相对运动，四路测距等于从地面站到中继卫星以及从中继卫星到用户卫星的距离之和。

TDRS 测距系统的测量噪声和系统误差通常小于 10m，这意味着用户航天器、中继卫星和地面站的运动必须建模到相似的精度。由于大多数由 TDRS 系统跟踪的用户卫星在几百公里的高度上绕地球运行，力模型至少应包含地球的谐波重力场以及大气阻力加速度。而地球同步中继卫星显著受月球和太阳的引力扰动以及太阳辐射压力的影响。为了避免不同卫星使用不同的模型，应考虑包含上述所有扰动的通用模型。建议使用到 20 阶的重力场模型来描述 500-1000 公里高度的用户卫星的运动，以达到所需的精度。

为了兼容地面站位置的建模，应考虑极移和真恒星时来描述地球自转。此外，在转换到通用惯性参考系时，还需考虑章动和岁差（例如 12000 年的平均赤道和春分点）。为了适当建模格林尼治时角，需要了解优于 0.01 秒的协调世界时（UT1）。

状态转移矩阵的变分方程可以基于仅包含重力场中二阶带谐波的简化力模型，与状态转移矩阵一起，描述状态向量相对于阻力和太阳辐射压力系数的偏导数的敏感度矩阵被集成到轨道确定中，以便调整这些参数。

2 TDRSOD 程序

TDRSOD 程序使用 TDRS 四路测距进行最小二乘法轨道确定。基于适当的测量集，可以同时调整单个用户航天器和最多两个中继卫星的轨道参数。跟踪数据提供在 TDRSOD.dat 文件中，其中包含每个时间的一条记录，结构如表 1 所述。单行头用来标记每一列，在读取时跳过。每行在测量的纪元后，给出地面站和航天器的识别号码以及四路中继测距。

表1 TDRSOD 跟踪数据文件的结构

补充的设置文件 TDRSOD.inp 提供了用户和 TDRS 卫星的先验状态向量和航天器相关参数以及相关的辅助信息。个别参数包括 TDRS 卫星的总数（≤ 2）和 ID、使用的地面站的总数（≤ 2）和 ID、所需的迭代次数、UT1-UTC 和 UTC-TAI 时间差、当前极点坐标，所有这些在输入文件的开头给出。在初始纪元之后，提供用户航天器和指定数量的 TDRS 卫星的状态向量及其先验标准差。航天器相关信息还包括每颗卫星的质量、面积、阻力系数和辐射压力系数，以及这些参数的先验不确定性。文件最后提供指定数量地面站的坐标。

每行提供一个从第 26 列开始的参数。输入时忽略初始字符，初始字符用于描述各个量的含义。以下列表给出了代表性的参数值：

TDRS 2 4 5
Stations 2 161 162
Iterations 4
UT1-UTC, UTC-TAI [s] +0.49 -32.00
x-pole,y-pole ["] -0.00651 +0.36588
Epoch (UTC) 1999/09/01 00:00:00.000
x UARS [m] 1476200.0 1000.0
y UARS [m] 5996200.0 1000.0
z UARS [m] -3209000.0 1000.0
vx UARS [m/s] -3854.0000 1.0
vy UARS [m/s] 3778.5000 1.0
vz UARS [m/s] 5302.2000 1.0
m [kg] , A [m²] 5968.3 27.22
CD, sigma(CD) 2.3 1.0
CR, sigma(CR) 1.3 0.1
x TDRS-4 [m] +20174293.6 1.0
y TDRS-4 [m] -37024903.8 1.0
z TDRS-4 [m] -982925.2 1.0
vx TDRS-4 [m/s] +2696.9634 0.001
vy TDRS-4 [m/s] +1471.5074 0.001
vz TDRS-4 [m/s] -100.5261 0.001
m [kg] , A [m²] 1668.9 40.0
CD, sigma(CD) 2.3 0.001
CR, sigma(CR) 1.3915 0.001
x TDRS-5 [m] -40783913.5 100.0
y TDRS-5 [m] 10622599.3 100.0
z TDRS-5 [m] 992633.1 100.0
vx TDRS-5 [m/s] -774.3896 0.1
vy TDRS-5 [m/s] -2976.0955 0.1
vz TDRS-5 [m/s] 18.8994 0.1
m [kg] , A [m²] 1718.4 40.0
CD, sigma(CD) 2.3 0.1
CR, sigma(CR) 1.4062 0.1
Sta WHSK/161 (xyz) [m] -1539385.74 -5160953.12 +3408202.16
Sta WH2K/162 (xyz) [m] -1539390.43 -5160968.83 +3408176.45

TDRSOD 程序使用 DE 多步法从指定的初始条件积分每个卫星的状态向量以及状态转移和敏感度矩阵。考虑到低地球用户卫星和地球同步中继卫星的轨道周期非常不同，分别独立积分各自的轨道，使用每个轨道最合适的积分步长。在光程时间迭代过程中使用多步插值。

最小二乘调整在指定的迭代次数内进行，每次迭代期间输出观测残差和计算的参数修正。为了保持简单的程序结构，不进行数据编辑或收敛性检查。所有估计参数需要提供先验标准差，应根据各状态向量分量或力模型参数的预期不确定性进行选择。（代码获取方式详见文末）

3 案例研究

在后续应用中，处理了 NASA 上层大气研究卫星（UARS）的 TDRS 测距数据，这些数据于 1999 年 9 月 1 日收集。数据集包括 14 批，每批 10-15 分钟的测量，每天均匀分布。其中，通过位于西经 41.0° 的 TDRS-4 卫星获得了三批数据，而其余的测量通过位于西经 174.3° 的 TDRS-5 卫星进行。所有数据都已预先校正了折射。

从地面跟踪独立得出的两颗 TDRS 卫星轨道作为 UARS 轨道确定的先验信息。相关状态向量和航天器参数如表 2所示。此外，表3提供了通过 TDRS-4 和 TDRS-5 进行四路测距测量的白沙地面站 WHSK 和 WH2K 天线的坐标。

为了适当进行最小二乘轨道确定，需要为状态向量分量以及阻力和太阳辐射压力系数指定先验标准差。在缺乏实际统计数据的情况下，可以通过考虑轨道特性、跟踪几何和数据分布获得适当的值。对于用户卫星，在大约 15 圈的轨道上具有充分覆盖的测量，预期可以从给定的测量中可靠确定其轨道元素。根据给定的数字，假定 UARS 卫星初始位置和速度的每个轴的不确定性分别为 1 公里和 1 米/秒。这些值不会对最终的最小二乘解构成严格约束，阻力系数的假定标准差为 1.0。另一方面，根据相关材料特性的认识，太阳辐射压力系数的先验标准差设定为 0.1。

对于地球同步高度的中继卫星，阻力不会引起任何轨道扰动，因此在轨道确定过程中无法校准相应的阻力系数。阻力系数的先验标准差仅用于避免正常方程的奇异性。实际上，先验值（2.3）在最小二乘调整中不会改变。对于其他参数，TDRS-4 跟踪数据仅限于少量数据，而 TDRS-5 跟踪数据覆盖全天。由于 TDRS-5 轨道的均匀采样，通过该卫星进行的四路测距测量可用于改进其轨道以及用户卫星的轨道。假设位置和速度的先验标准差分别为 100 米和 0.1 米/秒，太阳辐射压力系数的标准差为 0.1。这些数据考虑到 TDRS-5 轨道已经从独立跟踪数据中准确确定，同时避免用户卫星和 TDRS 卫星轨道参数之间的潜在相关性导致的不现实参数修正。最终，TDRS-4 卫星的先验标准差为 1 米（位置）、1 毫米/秒（速度）和 0.001（阻力系数和太阳辐射压力系数），本质上将其轨道限制在给定的先验轨迹。由于该卫星的 24 小时轨道周期数据覆盖不足，无法通过给定的 UARS 测距数据独立确定或改进其轨迹。

表2 UARS、TDRS-4 和 TDRS-5 在 1999/09/01 00:00 UTC 的先验轨道和航天器参数

表3 白沙天线位置

为了适当进行最小二乘轨道确定，需要为状态向量分量以及阻力和太阳辐射压力系数指定先验标准差。在缺乏实际统计数据的情况下，可以通过考虑轨道特性、跟踪几何和数据分布获得适当的值。对于用户卫星，在大约 15 圈的轨道上具有充分覆盖的测量，预期可以从给定的测量中可靠确定其轨道元素。根据给定的数字，假定 UARS 卫星初始位置和速度的每个轴的不确定性分别为 1 公里和 1 米/秒。这些值不会对最终的最小二乘解构成严格约束，阻力系数的假定标准差为 1.0。另一方面，根据相关材料特性的认识，太阳辐射压力系数的先验标准差设定为 0.1。

对于地球同步高度的中继卫星，阻力不会引起任何轨道扰动，因此在轨道确定过程中无法校准相应的阻力系数。阻力系数的先验标准差仅用于避免正常方程的奇异性。实际上，先验值（2.3）在最小二乘调整中不会改变。对于其他参数，TDRS-4 跟踪数据仅限于少量数据，而 TDRS-5 跟踪数据覆盖全天。由于 TDRS-5 轨道的均匀采样，通过该卫星进行的四路测距测量可用于改进其轨道以及用户卫星的轨道。假设位置和速度的先验标准差分别为 100 米和 0.1 米/秒，太阳辐射压力系数的标准差为 0.1。这些数据考虑到 TDRS-5 轨道已经从独立跟踪数据中准确确定，同时避免用户卫星和 TDRS 卫星轨道参数之间的潜在相关性导致的不现实参数修正。最终，TDRS-4 卫星的先验标准差为 1 米（位置）、1 毫米/秒（速度）和 0.001（阻力系数和太阳辐射压力系数），本质上将其轨道限制在给定的先验轨迹。由于该卫星的 24 小时轨道周期数据覆盖不足，无法通过给定的 UARS 测距数据独立确定或改进其轨迹。

表 4 1999/09/01 00:00 UTC UARS、TDRS-4 和 TDRS-5 的调整轨道和航天器参数

表4显示了 UARS、TDRS-4 和 TDRS-5 在 1999 年 9 月 1 日 00:00 UTC 的调整轨道和航天器参数的代表值，而相应的残差如图1所示。

总的来说，用户航天器的先验状态向量修正了大约 200 米和 0.1 米/秒，形式上的不确定性小了 5 到 10 倍。TDRS 卫星的初始状态向量基本保持不变，除了 TDRS-5 卫星的 z 分量被修改了 25 米。正如预期的那样，地球同步中继卫星的阻力系数完全没有受到影响，而用户卫星的太阳辐射压力系数在调整前后几乎相同。另一方面，用户卫星的阻力系数和 TDRS-5 的太阳辐射压力系数可以被良好地调整，并显著优于默认的先验值。

最终迭代中获得的残差显示了 5-10 米的总体测量和建模精度。显然，残差的分布不符合随机噪声的假设，而是表明存在系统误差。在缺乏用户卫星或 TDRS 卫星的独立跟踪数据的情况下，无法唯一地将这些误差归因于扰动力的不完整建模、预处理中介质修正的不完整考虑或测量过程中的系统误差。

4 具体使用说明（test_TDRSOD文件）

TDRSOD: 从跟踪和数据中继卫星测量中进行轨道确定

功能概述

TDRSOD（Tracking and Data Relay Satellite Orbit Determination）程序用于通过跟踪和数据中继卫星（TDRS）的测量数据进行轨道确定。该程序基于最小二乘法，处理卫星的四路测距数据来调整用户卫星和最多两个中继卫星的轨道参数。

主要参考文献

1. Montenbruck O., and Gill E., "Satellite Orbits: Models, Methods, and Applications," Springer Verlag, Heidelberg, Corrected 3rd Printing, 2005.
2. Vallado D. A., "Fundamentals of Astrodynamics and Applications," McGraw-Hill, New York, 4th edition, 2013.
3. 地球定向参数：http://celestrak.org/SpaceData/EOP-All.txt
4. 空间天气数据：http://celestrak.org/SpaceData/SW-All.txt
5. 行星历表：https://ssd.jpl.nasa.gov/planets/eph_export.html

使用说明

1. 准备工作

在运行test_TDRSOD.m脚本之前，确保以下文件在适当的目录中：

• DE440Coeff.mat：行星历表数据
• EOP-All.txt：地球定向参数数据
• SW-All.txt：空间天气数据

2. 设置全局变量

脚本设置了几个全局变量，用于存储天文常数、引力系数和其他辅助参数。

3. 加载所需数据

从DE440Coeff.mat文件中加载天文常数和行星系数。

SAT_Const
load DE440Coeff.mat
PC = DE440Coeff;

4. 读取地球引力场系数

从GGM03C.txt文件中读取地球引力场系数。

GGM03C_unnormalized

Cnm = zeros(361,361);
Snm = zeros(361,361);
fid = fopen('GGM03C.txt','r');
for n=0:360
    for m=0:n
        temp = fscanf(fid,'%d %d %f %f %f %f',[6 1]);
        Cnm(n+1,m+1) = temp(3);
        Snm(n+1,m+1) = temp(4);
    end
end
fclose(fid);

5. 读取地球定向参数

从EOP-All.txt文件中读取地球定向参数。

fid = fopen('EOP-All.txt','r');
while ~feof(fid)
    tline = fgetl(fid);
    k = strfind(tline,'NUM_OBSERVED_POINTS');
    if (k ~= 0)
        numrecsobs = str2num(tline(21:end));
        tline = fgetl(fid);
        for i=1:numrecsobs
            eopdata(:,i) = fscanf(fid,'%i %d %d %i %f %f %f %f %f %f %f %f %i',[13 1]);
        end
        for i=1:4
            tline = fgetl(fid);
        end
        numrecspred = str2num(tline(22:end));
        tline = fgetl(fid);
        for i=numrecsobs+1:numrecsobs+numrecspred
            eopdata(:,i) = fscanf(fid,'%i %d %d %i %f %f %f %f %f %f %f %f %i',[13 1]);
        end
        break
    end
end
fclose(fid);

6. 读取空间天气数据

从SW-All.txt文件中读取空间天气数据。

fid = fopen('SW-All.txt','r');
while ~feof(fid)
    tline = fgetl(fid);
    k = strfind(tline,'NUM_OBSERVED_POINTS');
    if (k ~= 0) 
        numrecsobs = str2num(tline(21:end));
        tline = fgetl(fid);
        for i=1:numrecsobs
            swdata(:,i) = fscanf(fid,'%i %d %d %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %f %i %i %f %i %f %f %f %f %f',[33 1]);
        end
        swdata = [swdata(1:27,:);swdata(29:33,:)];
        for i=1:4
            tline = fgetl(fid);
        end
        numrecspred = str2num(tline(28:end));
        tline = fgetl(fid);
        for i=numrecsobs+1:numrecsobs+numrecspred
            swdata(:,i) = fscanf(fid,'%i %d %d %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %f %i %i %f %f %f %f %f %f',[32 1]);
        end
        break
    end
end
fclose(fid);

7. 读取设置参数

从文件TDRSOD.txt中读取设置参数。

fid = fopen('TDRSOD.txt','r');

tline = fgetl(fid);
n_sat = str2num(tline(26:38));
n_sat = n_sat+1;

TdrsNo(2) = str2num(tline(39:51));
TdrsNo(3) = str2num(tline(52:64));

TdrsNo(1) = 0;

tline = fgetl(fid);
n_Sta = str2num(tline(26:38));

Sta(1).StaNo = str2num(tline(39:51));
Sta(2).StaNo = str2num(tline(52:64));

tline = fgetl(fid);
n_iterat = str2num(tline(26:38));

tline = fgetl(fid);
Y = str2num(tline(26:32));
M = str2num(tline(34:35));
D = str2num(tline(37:38));
h = str2num(tline(40:41));
m = str2num(tline(43:44));
s = str2num(tline(46:51));

Mjd0_UTC = Mjday(Y,M,D,h,m,s);
Aux(1).Mjd_UTC = Mjd0_UTC;
Aux(2).Mjd_UTC = Mjd0_UTC;
Aux(3).Mjd_UTC = Mjd0_UTC;

y0 = zeros(6,n_sat);
Sigy0 = zeros(6,n_sat);
p = zeros(2,n_sat);
Sigp = zeros(2,n_sat);

for i_sat=1:n_sat
    for i=1:6
        tline = fgetl(fid);
        y0(i,i_sat) = str2num(tline(26:38));
        Sigy0(i,i_sat) = str2num(tline(39:51));
    end
    
    tline = fgetl(fid);
    Aux(i_sat).mass = str2num(tline(26:38));
    Aux(i_sat).area_drag = str2num(tline(39:51));
    Aux(i_sat).area_solar = Aux(i_sat).area_drag;
    
    tline = fgetl(fid);
    Aux(i_sat).Cd = str2num(tline(26:38));
    Sigp(1,i_sat) = str2num(tline(39:51));
    
    tline = fgetl(fid);
    Aux(i_sat).Cr = str2num(tline(26:38));
    Sigp(2,i_sat) = str2num(tline(39:51));
    
    p(1,i_sat) = Aux(i_sat).Cd;
    p(2,i_sat) = Aux(i_sat).Cr;
end

Aux(1).n = 20;
Aux(1).m = 20;
Aux(1).n_a = 2;
Aux(1).m_a = 0;
Aux(1).n_G = 2;
Aux(1).m_G = 0;
Aux(1).sun = 1;
Aux(1).moon = 1;
Aux(1).sRad = 1;
Aux(1).drag = 1;
Aux(1).planets = 0;
Aux(1).SolidEarthTides = 0;
Aux(1).OceanTides = 0;
Aux(1).Relativity = 0;

Aux(2).n = 20;
Aux(2).m = 20;
Aux(2).n_a = 2;
Aux(2).m_a = 0;
Aux(2).n_G = 2;
Aux(2).m_G = 0;
Aux(2).sun = 1;
Aux(2).moon = 1;
Aux(2).sRad =

 1;
Aux(2).drag = 1;
Aux(2).planets = 0;
Aux(2).SolidEarthTides = 0;
Aux(2).OceanTides = 0;
Aux(2).Relativity = 0;

Aux(3).n = 20;
Aux(3).m = 20;
Aux(3).n_a = 2;
Aux(3).m_a = 0;
Aux(3).n_G = 2;
Aux(3).m_G = 0;
Aux(3).sun = 1;
Aux(3).moon = 1;
Aux(3).sRad = 1;
Aux(3).drag = 1;
Aux(3).planets = 0;
Aux(3).SolidEarthTides = 0;
Aux(3).OceanTides = 0;
Aux(3).Relativity = 0;

Rs = zeros(3,1);
for i_Sta=1:n_Sta
    tline = fgetl(fid);
    Rs(1) = str2num(tline(26:38));
    Rs(2) = str2num(tline(39:51));
    Rs(3) = str2num(tline(52:64));
    [Sta(i_Sta).lon,Sta(i_Sta).lat,Sta(i_Sta).h] = Geodetic(Rs);
end

fclose(fid);

8. 读取观测数据

从文件TDRSOD.dat中读取所有观测数据。

fid = fopen('TDRSOD.dat','r');
tline = fgetl(fid);
nobs = 0;
while ~feof(fid)
    nobs = nobs+1;
    tline = fgetl(fid);
    Y = str2num(tline(1:4));
    M = str2num(tline(6:7));
    D = str2num(tline(9:10));
    h = str2num(tline(13:14));
    m = str2num(tline(16:17));
    s = str2num(tline(19:24));
    Obs(nobs).Mjd_UTC = Mjday(Y,M,D,h,m,s);
    Obs(nobs).StaNo = str2num(tline(28:30));
    Obs(nobs).TdrsNo = str2num(tline(34));
    Obs(nobs).Range_4w = str2num(tline(37:46));
    Obs(nobs).Range_4w = Obs(nobs).Range_4w*1000;
end
fclose(fid);

9. 分配动态变量和对象

初始化轨道确定所需的动态变量和对象。

n_est = 8*n_sat;

X_apr = zeros(n_est,1);
dzdX = zeros(n_est,1);
dX = zeros(n_est,1);
SigX = zeros(n_est,1);

for i_sat=0:n_sat-1
    offset = 8*i_sat;
    for i=1:6
        X_apr(offset+i  ) = y0(i,i_sat+1);
    end
    for i=1:2
        X_apr(offset+6+i) = p(i,i_sat+1);
    end
    for i=1:6
        SigX(offset+i  )  = Sigy0(i,i_sat+1);
    end
    for i=1:2
        SigX(offset+6+i)  = Sigp(i,i_sat+1);
    end
end

P_apr = diag(SigX)*diag(SigX);

X = X_apr;

10. 运行迭代

通过迭代过程，逐步调整轨道参数，直至收敛。

for iterat=1:n_iterat
    R = zeros(n_est);
    d = zeros(n_est,1);

    LSQInit(n_est,X_apr-X,P_apr);
    n_obs = 1;
    Sum_res = 0;
    Sum_ressqr = 0;
    TrjData(1).y = y0(:,1);
    TrjData(2).y = y0(:,2);
    TrjData(3).y = y0(:,3);

    n_var = 8;
    TrjData(1).Y = zeros(7*n_var,1);
    TrjData(2).Y = zeros(7*n_var,1);
    TrjData(3).Y = zeros(7*n_var,1);

    for i=1:6
        TrjData(1).Y(i) = TrjData(1).y(i);
        TrjData(2).Y(i) = TrjData(2).y(i);
        TrjData(3).Y(i) = TrjData(3).y(i);
        for j=1:n_var
            if (i==j)
                TrjData(1).Y(6*j+i) = 1;
                TrjData(2).Y(6*j+i) = 1;
                TrjData(3).Y(6*j+i) = 1;
            else
                TrjData(1).Y(6*j+i) = 0;
                TrjData(2).Y(6*j+i) = 0;
                TrjData(3).Y(6*j+i) = 0;
            end
        end
    end

    for i=1:nobs
        for i_sat=1:n_sat
            AuxParam = Aux(i_sat);
            [TrjData(i_sat).y,TrjData(i_sat).Y] = Trj(Obs(i).Mjd_UTC,TrjData(i_sat).y,TrjData(i_sat).Y);
        end

        I_TDRS = -1;
        for i_sat=2:n_sat 
            if (TdrsNo(i_sat)==Obs(i).TdrsNo)
                I_TDRS = i_sat-1;
                break
            end
        end

        if (I_TDRS<0)
            continue
        end

        I_Sta = -1;
        for i_Sta=1:n_Sta
            if (Sta(i_Sta).StaNo==Obs(i).StaNo)
                I_Sta = i_Sta-1;
                break
            end
        end

        if (I_Sta<0)
            continue
        end

        [rho,drho_dx_User,drho_dx_TDRS] = Range_4W(Obs(i).Mjd_UTC,TrjData(1),TrjData(I_TDRS+1),Sta(I_Sta+1));

        res = Obs(i).Range_4w - rho;

        dzdX = zeros(n_est,1);
        offset = 8*I_TDRS;
        for k=1:8
            dzdX(k) = drho_dx_User(k);
            dzdX(offset+k) = drho_dx_TDRS(k);
        end

        n_obs = n_obs + 1;
        Sum_res = Sum_res + res;
        Sum_ressqr = Sum_ressqr + res*res;

        LSQAccumulate(n_est, dzdX, res, sigma_range);

        [year,mon,day,hr,min,sec] = invjday(Obs(i).Mjd_UTC);
        fprintf(' %4d/%2.2d/%2.2d  %2.2d:%2.2d:%6.3f', year, mon, day, hr, min, sec);
        fprintf('%6d', Obs(i).StaNo);
        fprintf('%6d', Obs(i).TdrsNo);
        fprintf('%13.4f', Obs(i).Range_4w/1000);
        fprintf('%12.4f', rho/1000);
        fprintf('%10.2f', Obs(i).Range_4w-rho);
        fprintf('\n');

        Aux(1).Mjd_UTC = Obs(i).Mjd_UTC;
        Aux(2).Mjd_UTC = Obs(i).Mjd_UTC;
        Aux(3).Mjd_UTC = Obs(i).Mjd_UTC;
    end
    Aux(1).Mjd_UTC = Mjd0_UTC;
    Aux(2).Mjd_UTC = Mjd0_UTC;
    Aux(3).Mjd_UTC = Mjd0_UTC;

    dX = LSQSolve(n_est, R, d);
    SigX = LSQStdDev(R, n_est);
    X = X + dX';

    for i_sat=0:n_sat-1
        offset = 8*i_sat;
        y0(:,i_sat+1) = X(offset+1:offset+6);
        Aux(i_sat+1).Cd = X(offset+7);
        Aux(i_sat+1).Cr = X(offset+8);
    end

    Mean = Sum_res/n_obs;
    RMS = sqrt(Sum_ressqr/n_obs-Mean*Mean);

    fprintf('%70s',' ');
    fprintf('_______\n');
    fprintf('%70s',' Mean ');
    fprintf('%.2f', Mean);
    fprintf(' m\n');
    fprintf('%70s',' RMS  ');
    fprintf('%.2f', RMS);
    fprintf(' m\n\n\n');
    fprintf('Results of iteration ');
    fprintf('%d', iterat);
   

 fprintf('\n\n');
    fprintf('  Parameter            ');
    fprintf('a priori   tot.corr. last corr.       final       sigma\n');

    for ii=1:8
        fprintf('  ');
        if (ii<4)
            fprintf('%s', Label(ii));
            fprintf(' [m]  ');
        elseif (3<ii && ii<7)
            fprintf('v');
            fprintf('%s', Label(ii-3));
            fprintf('[m/s]');
        elseif(6<ii && ii<9)
            fprintf('C');
            fprintf('%s', Label(ii-3));
            fprintf('     ');
        end
        fprintf(' User  ');
        fprintf('%14.2f', X_apr(ii));
        fprintf('%11.2f', X(ii)-X_apr(ii));
        fprintf('%11.2f', dX(ii));
        fprintf('%14.2f', X(ii));
        fprintf('%11.2f', SigX(ii));
        fprintf('\n');
    end

    for ii=9:16
        fprintf('  ');
        if (ii<12)
            fprintf('%s', Label(ii-8));
            fprintf(' [m]  ');
        elseif (11<ii && ii<15)
            fprintf('v');
            fprintf('%s', Label(ii-11));
            fprintf('[m/s]');
        elseif(14<ii && ii<17)
            fprintf('C');
            fprintf('%s', Label(ii-11));
            fprintf('     ');
        end
        fprintf(' TDRS-');
        fprintf('%d', TdrsNo(ceil(ii/8)));
        fprintf('%14.2f', X_apr(ii));
        fprintf('%11.2f', X(ii)-X_apr(ii));
        fprintf('%11.2f', dX(ii));
        fprintf('%14.2f', X(ii));
        fprintf('%11.2f', SigX(ii));
        fprintf('\n');
    end

    for ii=17:24
        fprintf('  ');
        if (ii<20)
            fprintf('%s', Label(ii-16));
            fprintf(' [m]  ');
        elseif (19<ii && ii<23)
            fprintf('v');
            fprintf('%s', Label(ii-19));
            fprintf('[m/s]');
        elseif(22<ii && ii<25)
            fprintf('C');
            fprintf('%s', Label(ii-19));
            fprintf('     ');
        end
        fprintf(' TDRS-');
        fprintf('%d', TdrsNo(ceil(ii/8)));
        fprintf('%14.2f', X_apr(ii));
        fprintf('%11.2f', X(ii)-X_apr(ii));
        fprintf('%11.2f', dX(ii));
        fprintf('%14.2f', X(ii));
        fprintf('%11.2f', SigX(ii));
        fprintf('\n');
    end
    fprintf('\n\n');
end

profile viewer
profile off

运行脚本时，MATLAB将自动执行这些步骤并输出结果。每次迭代后，会打印出轨道参数的调整情况和测量残差的统计数据。

Iteration 1

    Date         UTC        Sta   TDRS    obs [km]   comp [km]     o-c [m]
    1999/09/01  00:22: 1.000   162     5   79010.2586  79010.2461     12.46
     1999/09/01  00:22:31.000   162     5   78980.4438  78980.4297     14.09
     1999/09/01  00:23: 1.000   162     5   78954.4798  78954.4640     15.84
     1999/09/01  00:23:31.000   162     5   78932.3844  78932.3675     16.94
     1999/09/01  00:24: 1.000   162     5   78914.1732  78914.1544     18.79
     1999/09/01  00:24:31.000   162     5   78899.8554  78899.8346     20.84
     1999/09/01  00:25: 1.000   162     5   78889.4362  78889.4133     22.93
     1999/09/01  00:25:31.000   162     5   78882.9168  78882.8915     25.31
     1999/09/01  00:26: 1.000   162     5   78880.2936  78880.2657     27.85
     1999/09/01  00:26:31.000   162     5   78881.5586  78881.5282     30.39
     1999/09/01  00:27: 1.000   162     5   78886.6999  78886.6667     33.23
     1999/09/01  00:27:31.000   162     5   78895.7004  78895.6646     35.78
     1999/09/01  00:28: 1.000   162     5   78908.5399  78908.5013     38.62
     1999/09/01  00:28:31.000   162     5   78925.1932  78925.1517     41.54
     1999/09/01  00:29: 1.000   162     5   78945.6304  78945.5866     43.79
     1999/09/01  00:29:31.000   162     5   78969.8199  78969.7729     46.97
     ...
     1999/09/01  23:00:31.000   162     5   80196.9046  80183.8050  13099.58
     1999/09/01  23:01: 1.000   162     5   80327.1522  80313.8189  13333.26
     1999/09/01  23:01:31.000   162     5   80459.3313  80445.7790  13552.33
                                                                          _______
                                                                     Mean 5562.41 m
                                                                     RMS  5412.70 m
    
    Results of iteration 1
    
      Parameter            a priori   tot.corr. last corr.       final       sigma
      x [m]   User      1476200.00     -13.76     -13.76    1476186.24      12.95
      y [m]   User      5996200.00      62.72      62.72    5996262.72      11.51
      z [m]   User     -3209000.00     228.86     228.86   -3208771.14      17.30
      vx[m/s] User        -3854.00       0.02       0.02      -3853.98       0.01
      vy[m/s] User         3778.50      -0.12      -0.12       3778.38       0.02
      vz[m/s] User         5302.20       0.05       0.05       5302.25       0.01
      CD      User            2.30      -0.35      -0.35          1.95       0.10
      CR      User            1.30       0.00       0.00          1.30       0.10
      x [m]   TDRS-4   20174293.60       0.16       0.16   20174293.76       0.99
      y [m]   TDRS-4  -37024903.80      -0.10      -0.10  -37024903.90       1.00
      z [m]   TDRS-4    -982925.20      -0.00      -0.00    -982925.20       1.00
      vx[m/s] TDRS-4       2696.96       0.00       0.00       2696.96       0.00
      vy[m/s] TDRS-4       1471.51      -0.00      -0.00       1471.51       0.00
      vz[m/s] TDRS-4       -100.53      -0.00      -0.00       -100.53       0.00
      CD      TDRS-4          2.30       0.00      -0.00          2.30       0.00
      CR      TDRS-4          1.39      -0.00      -0.00          1.39       0.00
      x [m]   TDRS-5  -40783913.50      36.71      36.71  -40783876.79      13.56
      y [m]   TDRS-5   10622599.30      99.88      99.88   10622699.18      43.51
      z [m]   TDRS-5     992633.10      47.54      47.54     992680.64      55.16
      vx[m/s] TDRS-5       -774.39      -0.01      -0.01       -774.40       0.00
      vy[m/s] TDRS-5      -2976.10       0.00       0.00      -2976.09       0.00
      vz[m/s] TDRS-5         18.90      -0.01      -0.01         18.89       0.00
      CD      TDRS-5          2.30       0.00       0.00          2.30       0.10
      CR      TDRS-5          1.41      -0.01      -0.01          1.40       0.04