使用STK求出坐标转换矩阵

一、背景介绍

使用STK可以获得一颗卫星的地心惯性系J2000，地心固联坐标系，测站坐标系的x,y,z轴的指向分量，例如通过报告得到了测站坐标系（南东天）三个分量在J2000坐标系下的状态分量。如下图所示，那么从J2000坐标系到测站坐标系的坐标转换矩阵为

那么同样的通过测站坐标系在地心固联系下的状态分量，那么从地心固联系到测站坐标系的坐标转换矩阵为

那么根据坐标转换的传递性，可以求出J2000坐标系到地固系的坐标转换矩阵

二、数据实例

已知某卫星在历元时刻2024年3月17日 02：30：00-05：30：00每隔10s的位置和速度信息，state_Colum，根据STK数据求出的坐标转换矩阵，求出其轨道位置，与直接报表生成其位置进行对比。

（2.1）地固系到测站坐标系的转换

由于地固系和测站坐标系的关系，在115°东经，30°北纬的地方，从地固系到测站坐标系的转移矩阵为

Fixed2Topo=R_y(pi/2-30/180*pi)*R_z(115/180*pi)

求得的结果为

利用STK的报表生成，测站的南，东，天三个方向的坐标在地固系和J2000里的坐标。

clc;clear
%% 利用STK与matlab交互联合定轨
uiApplication = actxGetRunningServer('STK11.application');
% Get our IAgStkObjectRoot interface
root = uiApplication.Personality2;
checkempty = root.Children.Count;
if checkempty ~= 0
    root.CurrentScenario.Unload
    root.CloseScenario;
end
%% 根据你的需要设定场景的名称
root.NewScenario('OrbitDetermination');
StartTime = '17  Mar 2024 02:30:00.000';    % 场景开始时间
StopTime = '17 Mar 2024 05:30:00.000';     % 场景结束时间
root.ExecuteCommand(['SetAnalysisTimePeriod * "',StartTime,'" "',StopTime,'"']);
root.ExecuteCommand(' Animate * Reset');
root.ExecuteCommand(['New / */Facility Radar_01']);
root.ExecuteCommand(['SetPosition */Facility/Radar_01 Geodetic 30 115 0.0 MSL']);
Facility=root.CurrentScenario.Children.Item('Radar_01');
TopoCentric_Axes=Facility.VO.Vector.RefCrdns.Item(8);
TopoCentric_Axes.Visible=1;
TopoCentric_Axes.LabelVisible=1;
%% 记录每个时刻的测站坐标各个轴的在J2000坐标系下的分量(自己在STK新定义一个表格）
data3=root.ExecuteCommand(['Report_RM */Facility/Radar_01 Style "Vector_Topocentric" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep 10']);
%% 坐标变换矩阵
J20002Topop=zeros(3,3);
Sat_J2000=zeros(3,LineNum-2);
for i=1:LineNum-2
    %% 记录此时刻下的测站坐标对应的变换矩阵
    struct=regexp(data3.Item(i),',','split');
    J20002Topop(1,1)=str2num(struct{2});
    J20002Topop(1,2)=str2num(struct{3});
    J20002Topop(1,3)=str2num(struct{4});
    J20002Topop(2,1)=str2num(struct{5});
    J20002Topop(2,2)=str2num(struct{6});
    J20002Topop(2,3)=str2num(struct{7});
    J20002Topop(3,1)=str2num(struct{8});
    J20002Topop(3,2)=str2num(struct{9});
    J20002Topop(3,3)=str2num(struct{10});
    J2000Topo{i}=J20002Topop;
end
data4=root.ExecuteCommand(['Report_RM */Facility/Radar_01 Style "Vector_Topocentric_Fixed" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep 10']);
%% 坐标变换矩阵
Fixed=zeros(3,3);
Sat_Fixed=zeros(3,LineNum-2);
for i=1:LineNum-2
    %% 记录此时刻下的测站坐标对应的变换矩阵
    struct=regexp(data4.Item(i),',','split');
    Fixed(1,1)=str2num(struct{2});
    Fixed(1,2)=str2num(struct{3});
    Fixed(1,3)=str2num(struct{4});
    Fixed(2,1)=str2num(struct{5});
    Fixed(2,2)=str2num(struct{6});
    Fixed(2,3)=str2num(struct{7});
    Fixed(3,1)=str2num(struct{8});
    Fixed(3,2)=str2num(struct{9});
    Fixed(3,3)=str2num(struct{10});
    Fixed1{i}=Fixed;
end

其中Vector_Topocentric_Fixed的内容如图所示

Vector_Topocentric的内容如图所示

结果与通过旋转矩阵计算的结果一致，证明了三个坐标轴的分量按行排列即为从地心固联系到测站坐标系的转换矩阵。

（2.2）J2000坐标系到地固系的坐标转换

参考上面，我们同样采用两种方式生成坐标转换矩阵，第一个通过旋转坐标轴的方式

首先通过VGT模块生成J2000 X轴与Fixed X轴生成一个角度

通过插入一颗Planet，自己创建一个报告，内容即为theta角，求出02：30：00时刻的夹角为-147.65°，使用坐标转换矩阵，如下代码和结果如图所示

J20002Fixed=R_z(-147.65/180*pi)

采用计算测站坐标系在其上的分量的形式，求出测站坐标系的南，东，天三个方向在J2000坐标系的坐标分量为

得到从J2000坐标系到测站坐标系的坐标转换矩阵为

根据第一部分介绍的公式，从J2000到地心固联系的坐标转换即为

Fixed2To=[ 0.418981         ,        -0.269754          ,       -0.866999;
           0.539441         ,         0.842023          ,       -0.001296;
           0.730383         ,        -0.467151          ,        0.498308];
J20002To=[ -0.211309        ,          0.453154          ,       -0.866025;
                 -0.906308      ,           -0.422618        ,                0;
                 -0.365998      ,            0.784886        ,               0.5];
J20002Fixed=inv(Fixed2To)*J20002To

求出的结果为