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

ops/2024/11/2 0:35:46/

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文件）输入以下命令：

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

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

matlab">%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

matlab">%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的圆轨道，动力学模型为二体模型。
我们在分析问题时，绝大多数情况下都不会使用上述的默认参数，这就涉及到卫星轨道参数设置。卫星轨道参数设置，也存在多种方法，下面介绍：设置经典的轨道六根数。

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

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

轨道类型设置

matlab">kep.SizeShapeType = 'eSizeShapeAltitude';

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

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

+++

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

Orientation：轨道方向

matlab">kep.LocationType = 'eLocationTrueAnomaly';
kep.Orientation.AscNodeType = 'eAscNodeLAN';

+++

SizeShape：轨道的大小和形状

Orientation：轨道方向

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

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

matlab">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;

+++

显示卫星轨迹

matlab">sat.Propagator.InitialState.Representation.Assign(kep);
sat.Propagator.Propagate;
image-20220605185848825

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

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

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

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

image-20220605190632678

具体创建流程如下：

创建新的Facility并设置参数

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

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

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

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

创建新的Sensor并设置参数

创建新的Sensor的语句是：

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

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

image-20220605191856095

故代码如下：

matlab">senConstraints = sensor.AccessConstraints;altitude = senConstraints.AddConstraint('eCstrRange');%file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgSensor.html?Highlight=sensor.AccessConstraintsaltitude.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°

matlab">		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分析

matlab">        sat2rec=receiver.GetAccessToObject(Transmitter);sat2rec.ComputeAccess;end
end

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

<video id=“video” controls=""src=“https://pdswsj.oss-cn-beijing.aliyuncs.com/img/20220605_193050.mp4” preload=“none”>

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

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

使用STK逐一计算方法

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

image-20220605194445561

image-20220605194512589

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

image-20220605194719304

如上图所示。

使用MATLAB自动化获取report

正在学习，下周更新……

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

使用MATLAB自动化获取report

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

计算卫星可见性次数

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

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

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

获取卫星可见性数据

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

matlab">			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:lenacc.AER(lon*10+Lat).result(j+1+sum,:)={time(j) azimuth(j) range(j)};endsum=sum+len;end

实验结果

MATLAB结构体数据

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

image-20220608210756252

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

image-20220612181420914

实验效果视频

实验全流程视频如下

<video id=“video” controls=""src=“https://pdswsj.oss-cn-beijing.aliyuncs.com/img/20220612_181543.mp4” preload=“none”>

实验代码

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

matlab">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:9for lon=0:9if lon~=0 | Lat~=0facility = 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 terrainfacility.UseTerrain = true;% Set altitude to a distance above the groundfacility.HeightAboveGround = .05; % kmsensor = facility.Children.New('eSensor',['MySensor',num2str(lon),num2str(Lat)]);%         pattern1 = sensor.Pattern;%         sensor.FieldOfViewsenConstraints = sensor.AccessConstraints;altitude = senConstraints.AddConstraint('eCstrRange');%file:///E:/Program%20Files/AGI/STK%2011/Help/Programming/index.htm#DocX/STKObjects~IAgSensor.html?Highlight=sensor.AccessConstraintsaltitude.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:lenacc.AER(lon*10+Lat).result(j+1+sum,:)={time(j) azimuth(j) range(j)};endsum=sum+len;endendend
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

image-20220608210844209


http://www.ppmy.cn/ops/130275.html

相关文章

图书管理系统汇报

图书管理系统汇报

【1A536】图书管理系统汇报 项目介绍1.用户登录注册功能1. 1用户角色管理2.图书管理功能2.1 添加图书2.2 编辑图书2.3 删除图书 3.图书搜索和筛选3.1 图书搜索3.2 图书筛选 4.图书借阅、图书归还4.1 图书借阅4.2 图书归还 5.用户信息管理5.1上传头像5.2修改头像5.3 修改密码 项…
阅读更多...
虚拟机 Ubuntu 扩容

虚拟机 Ubuntu 扩容

文章目录 一、Vmware 重新分配 Ubuntu 空间二、Ubuntu 扩容分区 一、Vmware 重新分配 Ubuntu 空间 先打开 Vmware &#xff0c;选择要重新分配空间的虚拟机 点击 编辑虚拟机设置 &#xff0c;再点击 硬盘 &#xff0c;再点击 扩展 选择预计扩展的空间&#xff0c;然后点击 扩展…
阅读更多...
redis windows 7.0 下载

redis windows 7.0 下载

Redis 简介 Redis 是一个高性能的 key-value 数据库&#xff0c;广泛应用于缓存、消息队列、实时分析等场景。它支持多种数据结构&#xff0c;如字符串、哈希、列表、集合、有序集合等&#xff0c;并且提供了丰富的操作命令&#xff0c;能够满足各种复杂的数据处理需求。 下载…
阅读更多...
leetcode hot100【LeetCode 279. 完全平方数】java实现

leetcode hot100【LeetCode 279. 完全平方数】java实现

LeetCode 279. 完全平方数 题目描述 给定一个正整数 n&#xff0c;你需要找到若干个完全平方数&#xff08;比如 1, 4, 9, 16, ...&#xff09;&#xff0c;使得它们的和等于 n。你需要返回最少的完全平方数的个数。 示例 1: 输入: n 12 输出: 3 解释: 12 4 4 4示例 2…
阅读更多...
OB_GINS_day3

OB_GINS_day3

这里写目录标题 实现当前状态初始化实现预积分的初始化由于此时preintegration_options 是3&#xff08;也就是考虑odo以及earth rotation&#xff09;为预积分的容器添加需要积分的IMU积分因子接下来是添加新的IMU到preintegration中 实现当前状态初始化 这个state_curr的主要…
阅读更多...
[供应链] 邀请招标

[供应链] 邀请招标

1.邀请招标定义 邀请招标(Invitation to Bid by Request) 也称为有限竞争性招标(limited Competitive Bidding)或选择性招标(Selected Bidding) 邀请招标的采购方式下&#xff0c;采购人(如政府机构、企业或其他组织)不是公开发布招标信息&#xff0c;而是根据供应商或承包商…
阅读更多...
unity3d————三角函数练习题

unity3d————三角函数练习题

先上代码&#xff1a; public class SinCos : MonoBehaviour {public float moveSpeed 10f; //前进的速度public float changValue 5f; //左右的速度public float changeSize 5f; //左右的幅度float time 0;void Update(){this.transform.Translate(Vector3.forwa…
阅读更多...
java健身房管理系统源码(springboot)

java健身房管理系统源码(springboot)

项目简介 健身房管理系统实现了以下功能&#xff1a; 健身房管理系统实现了字典管理、健身房管理、教练管理、课程管理、器材管理、用户管理、管理员管理等功能。 &#x1f495;&#x1f495;作者&#xff1a;落落 &#x1f495;&#x1f495;个人简介&#xff1a;混迹java圈…
阅读更多...
最新文章

Copyright @ 2022~2023