MATLAB与STK互联：仿真并获取低轨卫星与指定区域地面站的可见性数据

MATLAB控制STK实现：仿真并获取低轨卫星与指定区域地面站的可见性数据

本次仿真主要参考了多篇文献和网站，包括但不限于：《Using MATLAB for STK Automation》、CSDN博文：拜火先知的博客_CSDN博客-笔记、AGI官网有关MATLAB的内容等内容，学习后很有收获，在此写下学习记录。

文章目录

MATLAB控制STK实现：仿真并获取低轨卫星与指定区域地面站的可见性数据
- 实验目的
- MATLAB与STK互联并创建场景
- 创建新卫星对象
- - 轨道类型设置
- 创建均匀分布的目标接收机
- - 创建新的Facility并设置参数
  - 创建新的Sensor并设置参数
  - 创建新的Receiver并设置参数
  - 创建Access分析
- 卫星相对于目标接收机的距离/速度/加速度/多普勒频率的计算
- - 使用STK逐一计算方法
  - 使用MATLAB自动化获取report
- 使用MATLAB自动化获取report
- - 计算卫星可见性次数
  - 获取卫星可见性数据
- 实验结果
- - MATLAB结构体数据
- 实验效果视频
- 实验代码

实验目的

使用MATLAB控制STK，实现如下功能：

生成1颗LEO卫星或产生一个LEO星座
在指定的区域内均匀分布目标接收机，指定区域大小为100km或500km
通过STK的access分析低轨卫星在过顶时，卫星相对于目标接收机的距离/速度/加速度/多普勒频率（分辨率为10s）

MATLAB与STK互联并创建场景

MATLAB与STK互联，主要有两种方式，一种是connect、一种是com口。这里主要采用com口形式，要比connect连接简单一些。
在MATLAB命令窗口（或新建.m文件）输入以下命令：

%这里改成自己STK的版本号
uiapplication = actxserver(‘STK11.application’);

键入下面的代码以创建具有自定义时间段的新场景。

%root是以后操作场景以及场景中对象的源头
root = app.Per	sonality2;
%新建场景
scenario = root.Children.New('eScenario','MATLAB_PredatorMission');
%将场景保存在指定目录下，同时可以将场景重新命名
root.SaveScenarioAs('E:\Program Files\STK\Documents\STK 11 (x64)\Scenario1_jz2');
scenario.SetTimePeriod('2 Jun 2022 16:00:00.000','3 Jun 2022 16:00:00.000');
scenario.StartTime = '2 Jun 2022 16:00:00.000';
scenario.StopTime = '3 Jun 2022 16:00:00.000';
root.ExecuteCommand('Animate * Reset');

创建新卫星对象

以下代码，将新建卫星对象。控制句柄为sat

%mysat是卫星在STK场景中的名字
sat = scenario.Children.New(‘eSatellite’,'mysat');
%执行下面命令，也可建立卫星。因为‘eSatellite’枚举序号为18
%sat = scenario.Children.New(18,‘mysat’);
%执行Propagate命令后，卫星将运行，即在STK中可以看到卫星运行轨迹
sat.Propagator.Propagate;

关于STK卫星对象属性值的设置，参见MATLAB与STK互联5：查看STK中对象的属性（2）—卫星对象属性梳理1_拜火先知的博客-CSDN博客

在新建卫星对象后，直接执行：sat.Propagate语句，会生成STK默认参数的卫星。默认参数为：倾角28.5°，轨道高度300km的圆轨道，动力学模型为二体模型。
我们在分析问题时，绝大多数情况下都不会使用上述的默认参数，这就涉及到卫星轨道参数设置。卫星轨道参数设置，也存在多种方法，下面介绍：设置经典的轨道六根数。

%将卫星轨道参数转换为开普勒六根数形式
kep=satPro.InitialState.Representation.ConvertTo(‘eOrbitStateClassical’);
%卫星轨道参数设置，可选择的形式有以下几种
%eOrbitStateClassical—经典轨道六根数形式
%eOrbitStateCartesian—位置、速度信息形式
%其他形式，大家可根据STK属性页查看，命名规则大致如上

下面是卫星参数设置，是本文重点

轨道类型设置

kep.SizeShapeType = 'eSizeShapeAltitude';

这里轨道形状类型，可以查看STK卫星轨道属性页面，有如下参数：

eSizeShapeSemimajorAxis—半长轴+偏心率型
eSizeShapeAltitude—远地点、近地点高度型
eSizeShapeRadius—远地点、近地点半径型
eSizeShapePeriod—轨道周期型
eSizeShapeMeanmotion—平动型

+++

LocationType：用于指定航天器在epoch轨道上的位置的元素

Orientation：轨道方向

kep.LocationType = 'eLocationTrueAnomaly';
kep.Orientation.AscNodeType = 'eAscNodeLAN';

+++

SizeShape：轨道的大小和形状

Orientation：轨道方向

Orientation.ArgOfPerigee：在卫星运动方向和轨道平面上，从上升节点到偏心矢量(轨道最低点)的角度。使用角维度。

Orientation.Inclination：轨道倾角。角动量矢量(垂直于轨道平面)和惯性Z轴之间的角。使用角维度。

kep.SizeShape.PerigeeAltitude = 500;
kep.SizeShape.ApogeeAltitude = 500;
kep.Orientation.Inclination = 60;
kep.Orientation.ArgOfPerigee = 0;
kep.Orientation.AscNode.Value = 0;
kep.Location.Value = 0;

+++

显示卫星轨迹

sat.Propagator.InitialState.Representation.Assign(kep);
sat.Propagator.Propagate;

通过matlab运行后，STK生成的“mysat”卫星如图。

创建均匀分布的目标接收机

在一个指定的目标区域内，按照每隔0.5个经度和0.5个纬度（每个约50km分布一个监测站），均匀分布着多个Facility、Sensor、Receiver（由于输出多普勒频率只能用transmitter到receiver，所以需要在Facility上附着Sensor和Receiver），这三者属于附着关系，Receiver附着在Sensor上，Sensor附着在Facility上。

首先附上生成的结果图，由图知，监测站等距离均匀地分布在一个正方形区域内，每隔约50km分布一个。

具体创建流程如下：

创建新的Facility并设置参数

首先需要进行for循环，在经纬度方向上个循环10次，也就是说一共在方圆500km内建立100个监测站。root.CurrentScenario.Children.New是建立一个新的Facility。

for Lat=0:9
    for lon=0:9
        facility = root.CurrentScenario.Children.New('eFacility', ['MyFacility',num2str(lon),num2str(Lat)]);

接着对Facility设置参数，值得一提的是，在对Facility进行位置定义时，facility.Position.AssignGeodetic每次增加0.5个经度/纬度。

    facility.Position.AssignGeodetic(41.9849+Lat*0.5,120.4039+lon*0.5,0) % Latitude, Longitude, Altitude
    % Set altitude to height of terrain
    facility.UseTerrain = true;
    % Set altitude to a distance above the ground
    facility.HeightAboveGround = .05; % km

创建新的Sensor并设置参数

创建新的Sensor的语句是：

sensor = facility.Children.New('eSensor',['MySensor',num2str(lon),num2str(Lat)]);

我们发现，如果单纯用该语句，那么sensor的视图会很难看，因此，我们需要关闭这个sight，但是由于找不到如何关闭，所以，只能接下来设置一下Constraints了，将Sensor的range_Max设置为0.

故代码如下：

senConstraints = sensor.AccessConstraints;
        altitude = senConstraints.AddConstraint('eCstrRange');%file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgSensor.html?Highlight=sensor.AccessConstraints
        altitude.EnableMax = true;
        altitude.Max = 0;
        % 可算出来了，设置了sensor的Constraints中的Range参数，让他为0，这样sensor就不会是个圆锥型了
        % IAgAccessConstraintCollection Collection
        % 对于sensor常规参数的修改，参考网址是：file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgAccessConstraintCollection.html?Highlight=eCstrLunarElevationAngle%20eCstrLighting
        %STK开发如果没有头绪，一个是查帮助，另一个就是用get、invoke查看相关信息。所以，可以上get

代码不多，但是很费力地才运行成功。

创建新的Receiver并设置参数

设置Receiver的最大接收范围是5000km，仰角为5°

		receiver = sensor.Children.New('eReceiver', ['MyReceiver',num2str(lon),num2str(Lat)]);
        %设置地面站约束
        facConstraints = receiver.AccessConstraints;
        %设置最大距离约束
        rangeCstr = facConstraints.AddConstraint('eCstrRange');
        rangeCstr.EnableMax = 1;
        rangeCstr.Max = 2000;
        %设置仰角约束
        elevationCstr = facConstraints.AddConstraint('eCstrElevationAngle');
        elevationCstr.EnableMin = 1;
        elevationCstr.Min = 5;

创建Access分析

        sat2rec=receiver.GetAccessToObject(Transmitter);
        sat2rec.ComputeAccess;
    end
end

最终生成的结果如下，视频中所展示的那样

卫星相对于目标接收机的距离/速度/加速度/多普勒频率的计算

获取的方法包括使用STK逐一计算和使用MATLAB自动化获取两种方法

使用STK逐一计算方法

选择：Analysis-Access，在Access for中选择“mysat”，任意选择一个MyFacility，这里拟选择“MyFacility00”，然后单击“Report & Graph Manager…”，在右边“Styles”中选择AER，可以生成详细报告和图表

此外，选择任意transmitter和receiver还可以产生多普勒频率的仿真值，使用的是“Link Budget - detailed”这个Styles

如上图所示。

使用MATLAB自动化获取report

正在学习，下周更新……

——————————————————————————————————————
更新：

使用MATLAB自动化获取report

实现MATLAB自动化需要以下几步：

计算卫星可见性次数

使用前文定义的sat2fac对象，根据STK提供的接口（IAgStkAccess Interface）：

 			acc_interval = sat2fac.ComputedAccessIntervalTimes;
            %可见次数如下
            count=acc_interval.Count;

计算出在一天时间内卫星可见次数。

获取卫星可见性数据

可见性数据包括各个时刻（分辨率为10秒），卫星相对于地面站的仰角，距离等要素。

			acc.AER(lon*10+Lat).result(1,:)={'Time' 'azimuth' 'range'};
            sum=0;
            for i=0:count-1
                [str,sto] = acc_interval.GetInterval(i);
                %结合上篇博文，获取第一个可见弧段内的方位角
                aerDP = sat2fac.DataProviders.Item('AER Data').Group.Item('VVLH CBF').Exec(str,sto,10);
                azimuth = cell2mat(aerDP.DataSets.GetDataSetByName('Azimuth').GetValues);
                range = cell2mat(aerDP.DataSets.GetDataSetByName('Range').GetValues);
                time = aerDP.DataSets.GetDataSetByName('Time').GetValues;
                len=length(range);
                for j=1:len
                    acc.AER(lon*10+Lat).result(j+1+sum,:)={time(j) azimuth(j) range(j)};
                end
                sum=sum+len;
            end

实验结果

MATLAB结构体数据

我们讲个地面站与卫星的参数保存为一个结构体，字段1代表Facility01地面站相对于卫星Sat的各个可见时刻的数据信息。如图所示：

该图就是Facility01地面站相对于卫星Sat的各个可见时刻的具体数据信息。

实验效果视频

实验全流程视频如下

实验代码

写了一晚上程序，可以实现这个目标了，先附上最终代码

clear all
clc
app = actxserver('STK11.application');
root = app.Personality2;
scenario = root.Children.New('eScenario','MATLAB_PredatorMission');
root.SaveScenarioAs('E:\Program Files\STK\Documents\STK 11 (x64)\Scenario1_jz2');
scenario.SetTimePeriod('2 Jun 2022 16:00:00.000','3 Jun 2022 16:00:00.000');
scenario.StartTime = '2 Jun 2022 16:00:00.000';
scenario.StopTime = '3 Jun 2022 16:00:00.000';
root.ExecuteCommand('Animate * Reset');
%facility = scenario.Children.New('eFacility','GroundStation'); 
% areaTarget = root.CurrentScenario.Children.New('eAreaTarget', 'MyAreaTarget');
% root.BeginUpdate();
% areaTarget.AreaType = 'eEllipse';
% ellipse = areaTarget.AreaTypeData;
% %ellipse.Add(48.897, 18.637);
% ellipse.SemiMajorAxis = 85.25; % in km (distance dimension)
% ellipse.SemiMinorAxis = 85.25; % in km (distance dimension)
% ellipse.Bearing = 44; % in deg (angle dimension)
% root.EndUpdate();
sat = scenario.Children.New(18,'mysat');
kep = sat.Propagator.InitialState.Representation.ConvertTo('eOrbitStateClassical');
kep.SizeShapeType = 'eSizeShapeAltitude';
kep.LocationType = 'eLocationTrueAnomaly';
kep.Orientation.AscNodeType = 'eAscNodeLAN';
kep.SizeShape.PerigeeAltitude = 500;
kep.SizeShape.ApogeeAltitude = 500;
kep.Orientation.Inclination = 60;
kep.Orientation.ArgOfPerigee = 0;
kep.Orientation.AscNode.Value = 0;
kep.Location.Value = 0;
sat.Propagator.InitialState.Representation.Assign(kep);
sat.Propagator.Propagate;
sensor = sat.Children.New('eSensor','MySensor_sat');
Transmitter = sensor.Children.New('eTransmitter', 'Mytransmitter');
Transmitter.SetModel('GPS Satellite Transmitter Model');
%建立卫星对象，并设置参数
% satellite = root.CurrentScenario.Children.New('eSatellite', 'MySatellite');
% scenPath = 'E:\Program Files\STK\Documents\STK 11 (x64)\Scenario1_learn3';
% satelliteFilePath = [scenPath '\ISS_25544.sa'];
% satellite.ExportTools.GetEphemerisStkExportTool.Export(satelliteFilePath);

% IAgFacility facility: Facility Object
for Lat=0:9
    for lon=0:9
        if lon~=0 | Lat~=0
            facility = root.CurrentScenario.Children.New('eFacility', ['MyFacility',num2str(lon),num2str(Lat)]);
            facility.Position.AssignGeodetic(41.9849+Lat*0.5,120.4039+lon*0.5,0) % Latitude, Longitude, Altitude
            % Set altitude to height of terrain
            facility.UseTerrain = true;
            % Set altitude to a distance above the ground
            facility.HeightAboveGround = .05; % km
            sensor = facility.Children.New('eSensor',['MySensor',num2str(lon),num2str(Lat)]);
            %         pattern1 = sensor.Pattern;
            %         sensor.FieldOfView
            senConstraints = sensor.AccessConstraints;
            altitude = senConstraints.AddConstraint('eCstrRange');%file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgSensor.html?Highlight=sensor.AccessConstraints
            altitude.EnableMax = true;
            altitude.Max = 0;
            % 可算出来了，设置了sensor的Constraints中的Range参数，让他为0，这样sensor就不会是个圆锥型了
            % IAgAccessConstraintCollection Collection
            % 对于sensor常规参数的修改，参考网址是：file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgAccessConstraintCollection.html?Highlight=eCstrLunarElevationAngle%20eCstrLighting
            %STK开发如果没有头绪，一个是查帮助，另一个就是用get、invoke查看相
            %关信息。所以，上get
            %         pattern1.ConeAngle = 85;
            receiver = sensor.Children.New('eReceiver', ['MyReceiver',num2str(lon),num2str(Lat)]);
            %设置地面站约束
            facConstraints = receiver.AccessConstraints;
            %设置最大距离约束
            rangeCstr = facConstraints.AddConstraint('eCstrRange');
            rangeCstr.EnableMax = 1;
            rangeCstr.Max = 2000;
            %设置仰角约束
            elevationCstr = facConstraints.AddConstraint('eCstrElevationAngle');
            elevationCstr.EnableMin = 1;
            elevationCstr.Min = 5;
            sat2fac=sat.GetAccessToObject(facility);
            sat2fac.ComputeAccess;
            %         aerDP = sat2fac.DataProviders.Item('AER Data').Group.Item('VVLH CBF').Exec(scenario.StartTime,scenario.StopTime,60);
            %         azimuth = cell2mat(aerDP.DataSets.GetDataSetByName('Range').GetValues);
            acc_interval = sat2fac.ComputedAccessIntervalTimes;
            %可见次数如下
            count=acc_interval.Count;
            %可见的弧段起始时间列表
            %acc_interval.ToArray(0,-1)
            %str、sto分别是第一个弧段的起始时间，是字符串类型。
            acc.AER(lon*10+Lat).result(1,:)={'Time' 'azimuth' 'range'};
            sum=0;
            for i=0:count-1
                [str,sto] = acc_interval.GetInterval(i);
                %结合上篇博文，获取第一个可见弧段内的方位角
                aerDP = sat2fac.DataProviders.Item('AER Data').Group.Item('VVLH CBF').Exec(str,sto,10);
                azimuth = cell2mat(aerDP.DataSets.GetDataSetByName('Azimuth').GetValues);
                range = cell2mat(aerDP.DataSets.GetDataSetByName('Range').GetValues);
                time = aerDP.DataSets.GetDataSetByName('Time').GetValues;
                len=length(range);
                for j=1:len
                    acc.AER(lon*10+Lat).result(j+1+sum,:)={time(j) azimuth(j) range(j)};
                end
                sum=sum+len;
            end
        end
    end
end

% satellite = root.GetObjectFromPath('/Satellite/mysat');
% receiver = root.GetObjectFromPath('/Facility/MyFacility99');
% senConstraints = sensor.AccessConstraints;
% altitude = senConstraints.AddConstraint('eCstrAltitude');%file:///E:/Program%20
% % Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgSensor.html?Highlight=sensor.AccessConstraints
% altitude.EnableMax = true; 
% altitude.Max = 2;
% AER()
% 对于sensor常规参数的修改，参考网址是：file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgAccessConstraintCollection.html?Highlight=eCstrLunarElevationAngle%20eCstrLighting

各个站点的结果输出在了一个结构体中，每个结构体包含以下三种元素（time azimuth range），其中range就是星地间距离，可以作为整个系统的输出。数据保存在这个链接，下载即可查看https://pdswsj.oss-cn-beijing.aliyuncs.com/img/matlab.mat