一文教你做MATLAB与STK编程（MATLAB与STK混合编程实例，含代码）

简介：通过STK的连接模块（connectors），实现STK与应用程序Matlab的交互。单独一方干不成或难度大，通过互联可充分发挥STK与第三方应用程序各自的优势。STK与Matlab互联是两个强大成熟商业软件的强强联合。STK提供超过200个Matlab格式化过封装好的命令。这些接口函数能够大大提升编程效率。

STK强项：轨道计算功能丰富、 可靠。
STK弱项：STK自身无法通过编程实现对某些复杂航天任务的仿真分析。例如：循环计算、 复杂的嵌套迭代、 复杂的收敛判据。

Matlab强项：使用便捷、 能编程实现复杂逻辑。

Matlab弱项：没有天生的轨道计算能力。

一、本博客的MATLAB与STK互联实例任务：

背景：给定一个四颗卫星组成的星座，计算分析两天时间内该星座对全球每一点的最大重访时间。

任务：计算星座对地面点的最大重访时间。

二、直接上结果图：

1.  该动图为MATLAB与STK互联的直观显示，MATLAB运行程序后，与STK建立connectors，自动打开STK工程，并自动建立对应场景，通过遍历地球表面地面站，计算重访时间，循环得到两天时间内该星座对全球每一点的最大重访时间。

2.  下图为MATLAB显示的waitbar，可以自动显示循环进度以及预计剩余时间。

3. 下图为四星星座对全球每一点的最大重访时间二维图以及三维图。纬度变化范围为：-90°~90°；经度的变化范围为：0°~360°。

步进角度为5度的结果：

步进角度为2度的结果：  

三、MATLAB与STK互联步骤

MATLAB与STK互联主要可以通过电脑串口COM或者STK的connectors链接器进行互联。本文主要介绍connection链接器方法。该部分已经有很多博主进行分享，本博客进行简要介绍。建议先安装MATLAB2018b，再安装STK11.6。这两个版本的软件能够自动建立connectors，完成软件的互联。

1. 软件版本截图：

2. 链接connection步骤：（注意：前提是已经安装好MATLAB2018b） 

① 以管理员身份打开STK11.6的install.exe。界面如下：

红色框内选项全部勾选，由于我电脑已经安装好了STK11.6，因此我的红色框内选项是灰色的，没有安装STK11.6的电脑是可以勾选的。

② 安装好STK11.6后，打开STK，查看STK和Matlab互联连接器connectors的状态。点击STK菜单栏的edit选项卡，继续点击Preferences，可以看到与STK连接的MATLAB版本，下图代表STK与MATLAB互联完成。

③ 然后再打开MATLAB，输入stkInit（区分大小写），对STK的连接进行初始化，可以看到如下提示，这就表示MATLAB已经与STK进行连接，此后就可以利用MATLAB对STK进行编程。

 3. 航天应用常用接口函数的使用。

① 初始化需要的接口函数

stkInit—— 完成STK和Matlab的互联

conid=stkOpen(stkDefaultHost); —— 返回互联成功的主机端口的连接句柄

② 初始窗口管理

if  stkValidScen == 1

stkUnload('/*')

end —— 如果已经有打开的场景，则关闭场景。

③ 建立场景

stkNewObj('/','Scenario','场景名称'); —— 建立给定名称的场景。

stkSetTimePeriod('10 Apr 2024 00:00:00.0','12 Apr 2024 00:00:00.0','GREGUTC');

—— 设置场景的起止时间和采用的时间系统。

stkSetEpoch('10 Apr 2024 00:00:00.0','GREGUTC'); —— 设置场景的历元。

stkSyncEpoch; —— 同步aeroToolbox和STK场景历元。

rtn = stkConnect(conid,'Animate','Scenario/场景名称','SetValues "10 Apr 2024

00:00:00.0" 60 0.1'); —— 设置STK场景动画历元。

rtn = stkConnect(conid,‘Animate’,‘Scenario/场景名称’,‘Reset’); —— 设置动画时

间复位。

④ 建立航天器

stkNewObj('*/','Satellite','航天器名称'); —— 建立卫星。

stkSetPropClassical('objPath', 'propagator', 'coordSystem', tStart, tStop, dt,

orbitEpoch, semimajorAxis, eccentricity,inclination, argOfPerigee, RAAN,

meanAnomaly, coordEpoch)

⑤ 关闭连接

stkClose(conid);

4. 所有STK与MATLAB互联语法如下PDF显示：

四、MATLAB与STK互联主要代码如下，如有代码问题，请加V（Lucky_YQW）交流。

clc
clear all
close all
stkClose('ALL')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算星座对地面点的最大重访时间
%详细背景：给定一个星座，计算分析两天时间内该星座对全球每一点的最大重访时间。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 互联
stkInit; % 建立连接
conid=stkOpen; % 得到连接句柄（用于发送指令）
delete(get(0,'children')); % 关闭其他绘图窗口
scen_open = stkValidScen; % 如果有已打开的场景，则自动关闭
if scen_open == 1
    stkUnload('/*');
end

%% 建立场景
% 新建场景
stkNewObj('/','Scenario','Test');
% 设定时间
stkSetTimePeriod('26 Apr 2024 00:00:00.00','28 Apr 2024 00:00:00.00','GREGUTC');
% 使aeroToolbox中函数的历元时间与场景一致
stkSetEpoch('26 Apr 2024 00:00:00.00','GREGUTC');
stkSyncEpoch;

%% 新建四颗卫星
stkNewObj('*/','Satellite','sat1');
stkNewObj('*/','Satellite','sat2');
stkNewObj('*/','Satellite','sat3');
stkNewObj('*/','Satellite','sat4');
% 设置起止时间、历元时刻、积分步长
t_start=0;
t_stop=1*(24*3600);   % 两天的时间
dt=60;
orbitEpoch=t_start;
% 基准星轨道根数
a_MB=7000*1000;
e_MB=0;
i_MB=60*pi/180;
w_MB=0*pi/180;
Raan_MB=0*pi/180;
M_MB=0*pi/180;

% 输入并积分星座轨道
stkSetPropClassical('*/Satellite/sat1','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_MB,e_MB,i_MB,w_MB,Raan_MB,M_MB);
stkPropagate('*/Satellite/sat1', t_start, t_stop)
stkSetPropClassical('*/Satellite/sat2','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_MB,e_MB,i_MB,w_MB,Raan_MB+pi/4,M_MB+pi);
stkPropagate('*/Satellite/sat2', t_start, t_stop)
stkSetPropClassical('*/Satellite/sat3','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_MB,e_MB,i_MB,w_MB,Raan_MB+pi/2,M_MB+2*pi);
stkPropagate('*/Satellite/sat3', t_start, t_stop)
stkSetPropClassical('*/Satellite/sat4','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_MB,e_MB,i_MB,w_MB,Raan_MB+3*pi/4,M_MB+pi/2);
stkPropagate('*/Satellite/sat4', t_start, t_stop);


%% 建立地面站，循环计算不同位置的最大重访时间
% 新建地面站
stkNewObj('*/','Facility','Point');    % Phi; Lamda 分别是纬度，经度，纬度变化范围是-90到90度，经度0-360度
% 地面站参数
% 经度参数
Lamda_min = 0;
delta_Lamda = 5;
n_Lamda = 360/delta_Lamda;
% 纬度参数
Phi_min = -90;
delta_Phi = 5;
n_Phi = 180/delta_Phi;
% 循环计算
h=waitbar(0,'开始地面站扫描');
for i=1:n_Lamda
    for j=1:n_Phi
        tic
        % 设置地面站经纬度坐标
        namda = (Lamda_min+(i-1)*delta_Lamda)/180*pi;
        X(i,j)=namda/pi*180;
        phi = (Phi_min+(j-1)*delta_Phi)/180*pi;
        Y(i,j)=phi/pi*180;
        stkSetFacPosLLA('Scenario/Test/Facility/Point',[phi; namda; 0]);
        % 根据对地观测需求设置地面点最低跟踪仰角
        stkConnect(conid,'SetConstraint','Scenario/Test/Facility/Point','ElevationAngle Min 50');
        %计算各星访问时间

        %%%%%%%%%%%%%此处省略%%%%%%%%%%%
        
        % 计算最大重访时间，直至循环结束
        temp=0;
        if size(interval1)==[0 0] & size(interval2)==[0 0] & size(interval3)==[0 0] & size(interval4)==[0 0]
            Z(i,j)=0;
        else
            if size(interval1)~=[0 0]
                temp=[temp,interval1.start,interval1.stop];
            end
            if size(interval2)~=[0 0]
                temp=[temp,interval2.start,interval2.stop];
            end
            if size(interval3)~=[0 0]
                temp=[temp,interval3.start,interval3.stop];
            end
            if size(interval4)~=[0 0]
                temp=[temp,interval4.start,interval4.stop];
            end

        %%%%%%%%%%%%%此处省略%%%%%%%%%%%

        end
        toc
        flag_num = (i-1) * n_Phi + j;
        if i == 1 && j == 1
        per_time = toc;
        end
        time_remain = per_time * (n_Lamda*n_Phi - flag_num);
        str=['扫描进度：',num2str(100*flag_num/n_Lamda/n_Phi),'%','    预估剩余时间:',num2str(time_remain),'秒'];       
        waitbar(flag_num/n_Lamda/n_Phi,h,str);
    end
end
delete(h);  %close(h)
save('res_inteval_5.mat','X','Y','Z');

figure
surf(X,Y,Z);
xlabel('经度（度）');
ylabel('纬度（度）');
zlabel('最大重访时间间隔（小时）');
colorbar
xlim([0 355]);
ylim([-90 85]);

%% 关闭连接
stkClose(conid);