一次Python与STK12.2联合仿真心得体会

        本人在做一个项目的时候,需要使用STK12.2创建大量卫星,并给每个卫星设置成像载荷,计算整个场景时间下,所有成像载荷与地面目标的Access数据,最后还要将数据导出用于其他处理计算。

        在尝试手动使用STK12.2实现以上效果失败后,考虑使用Python对STK12.2编程实现。本人编程语言方面只学过Python语言(而且学得很普通),整个联合仿真过程下来,感觉就是:入门较难、入门后感觉较简单;二次开发的水平取决于对STK12.2这个软件了解的多少,了解的越多,编程实现某些功能就越容易。

        本次联合仿真大致过程:

(一)软件准备:

        STK12.2是在某宝上花钱买的。我个人在安装软件上,更偏向于能用钱解决的就用钱解决,无论是商家远程安装还是自己按照商家的步骤安装,效率都更高,而自己从网上找免费的渠道安装软件费时费力还不一定能成功。Python是自己按照版本对应关系下载的,我使用的Python版本是3.10.9。我是在Pycharm Community Edition 2023.1.1上进行编程。

(二)STK12与Python开发的环境配置:

        互联教程是在CSDN上搜到的是作者“一只大佬铁”的教程,链接我放在下面了。

STK12与Python联合仿真(一):环境搭建-CSDN博客

        如果按照他的步骤完整走下来,最后可以在STK中打开一个Jupyter notebook进行编程。我按照步骤来搞,没有完全成功,但是Pycharm已经可以对STK12.2进行编程了,具体什么原因我也不太清楚。环境配置这块我搞得也是稀里糊涂,建议大家电脑上安装一个Pycharm。然后按照这位作者的步骤试一试。

STK12.2+Python开发(一):STK和Python开发环境配置_python链接stk-CSDN博客

        如果已经成功连接了,那么STK12.2自带的例子程序是可以直接运行的。正常情况下,例子程序存在路径C:\Program Files\AGI\STK 12\CodeSamples\Automation和路径C:\Program Files\AGI\STK 12\CodeSamples\CustomApplications。这两个文件夹中不仅存有python例子程序,还有C++、Java等程序。

(三)编程应用        

        安装与互联完成后,更头疼的就来了。初次进行互联编程真的就感觉无从下手,没有系统的教程,能查到的一些文章页很少,而且很多还是STK11版本的。除了CSDN上可以搜到的教程,这里也是推荐知乎作者“小脏鱼儿”的教程,虽然教程是STK11版本的,但是其中关于对象模型和模型调用的讲解非常值得学习,也是适用于STK12的,链接放下面了。

Python与STK联合仿真-第三讲 - 知乎 (zhihu.com)

        除了例子程序、网上的教程以外,能够参考的资料就是STK自带的非常完整的英文说明书,Programming Interface Help,在STK中可以直接打开。

        打开后可以在里面直接搜索对象和接口的使用方法,但是首先要知道有什么对象、有什么接口,他们之间的关系是什么。这就需要在说明书中下载另一个关键的资料——接口关系图。

        可以在说明书中看到,STK Object Model的库有STK Objects、STK X、STK VGT等6个。每个库都有如下的一张接口关系图,里面包含了该库所有的接口、对象、对象的属性和方法以及他们之间的关系。最常用的就是STK Objects的关系图,这个关系图满足的了我要实现的所有功能。

        有了说明书和关系图,再加上例子程序和网上能搜到的教程,通过自己的慢慢摸索,就可以一步一步掌握STK12与Python二次开发的方法。目前我总结了以下开发流程:

第一步:明确自己要使用python与STK12.联合仿真达成什么效果

第二步:明确哪些效果是使用STK12实现的,哪些是python实现的。以构建2000颗卫星为例,那么STK12要实现构建一颗卫星的效果,python实现将这个过程循环2000次。

第三步:将第二步STK12要实现的效果,先在该软件中实现。如果在软件中都不知道如何实现,那么也很难编程实现。例如,先在STK软件中构建一颗卫星。

第四步:查找关系图,明确软件中操作STK实现的功能需要使用那些接口、对象。例如构建一颗卫星需要IAgStkObjectRoot、IAgStkRoot、IAgSatellite。再使用说明书,搜索这些接口,查看接口下的方法,对象以及属性,哪些可以实现这些功能,同时也要配合关系图,查看是否符合接口之间的逻辑关系,如果都符合,即可以使用该对象以及方法。可以点击这些对象、方法,它们的具体用法会有显示。

下图是IAgSTKObject接口包含的方法和属性。

下图是该接口中Children在各种语言中的具体用法

        如果报错了,很有可能就是接口找错了,或者接口之间关系搞错了。总之找到可以实现功能的接口是关键,从结构图上找一个满足功能的接口是一个漫长的过程;如果接口找对了,但是接口之间关系搞错了,大概率是因为不熟悉结构之间关系,多试几次就熟练了。

终极秘笈:

        实在不行就使用AI写代码,这个我已经尝试过,效果很好。只需要告诉AI,使用Python语言编程,使用STK12达成的效果即可。如果有报错,那么就把报错语句发给AI,他会根据报错进一步修改代码。

下面这段代码是根据我提出的要求(输入卫星TLE数据,可以得到相应的轨道六要素。),AI给我的第一版代码,在我的电脑上运行是有问题的。

from skyfield.api import EarthSatellite, load
from math import degrees

def tle_to_orbital_elements(tle_line1, tle_line2):
    # 创建卫星对象
    satellite = EarthSatellite(tle_line1, tle_line2)
    
    # 获取轨道要素
    semi_major_axis = satellite.model.a  # 半长轴 (km)
    eccentricity = satellite.model.e  # 偏心率
    inclination = degrees(satellite.model.inclination)  # 倾角 (度)
    raan = degrees(satellite.model.raan)  # 升交点赤经 (度)
    arg_of_perigee = degrees(satellite.model.argument_of_perigee)  # 近地点幅角 (度)
    mean_anomaly = degrees(satellite.model.mean_anomaly)  # 平近点角 (度)
    
    return {
        "Semi-major axis (km)": semi_major_axis,
        "Eccentricity": eccentricity,
        "Inclination (deg)": inclination,
        "RAAN (deg)": raan,
        "Argument of Perigee (deg)": arg_of_perigee,
        "Mean Anomaly (deg)": mean_anomaly
    }

# 主程序
if __name__ == "__main__":
    print("请输入TLE数据:")
    tle_line1 = input("第一行: ")
    tle_line2 = input("第二行: ")
    
    orbital_elements = tle_to_orbital_elements(tle_line1, tle_line2)
    
    print("\n轨道六要素:")
    for key, value in orbital_elements.items():
        print(f"{key}: {value:.6f}")

        经过多次报错反馈以及提出修改建议,AI最终给我的代码如下。AI知道要用什么库,但是他可能对库中的具体结构不清楚,通过报错反馈,他慢慢地推测出了库的结构,给出了最有可能正确的程序。

from skyfield.api import EarthSatellite
import math

# 假设你已经创建了一个 EarthSatellite 对象
satellite = EarthSatellite(line1, line2, name="Example Satellite")

# 获取轨道参数
model = satellite.model

# 访问轨道参数
inclination = math.degrees(model.inclo)
eccentricity = model.ecco
right_ascension = math.degrees(model.nodeo)
argument_of_perigee = math.degrees(model.argpo)
mean_motion = model.no  # 注意:这可能是弧度/分钟
mean_anomaly = math.degrees(model.mo)
semi_major_axis = model.a  # 添加半长轴

# 转换平均运动为每天圈数
mean_motion_revs_per_day = mean_motion * 1440 / (2 * math.pi)

print(f"Inclination: {inclination:.2f} degrees")
print(f"Eccentricity: {eccentricity:.6f}")
print(f"Right Ascension of Ascending Node: {right_ascension:.2f} degrees")
print(f"Argument of Perigee: {argument_of_perigee:.2f} degrees")
print(f"Mean Motion: {mean_motion_revs_per_day:.6f} revolutions per day")
print(f"Mean Anomaly: {mean_anomaly:.2f} degrees")
print(f"Semi-major Axis: {semi_major_axis:.2f} km")  # 添加半长轴的输出

# 额外信息
print(f"Epoch Year: {model.epochyr}")
print(f"Epoch Days: {model.epochdays:.8f}")
print(f"B* drag term: {model.bstar:.8f}")

(四)例子程序

        下面给出我例子程序,实现功能:读取TLE数据构建卫星,以及读取传感器信息文件(该文件为我自己设定的)为每颗卫星设置传感器。这段程序已经可以实现文章开头所说的一部分功能了,至于完整的功能(构建地面目标,计算Access数据),我准备放到下一篇文章中详细说明。

import os
from skyfield.api import EarthSatellite
import math
import openpyxl
#######
#这段程序必须要有,算是一个初始化程序,可以在每个里程序中看到。
try:
    if os.name == "nt":
        from agi.stk12.stkdesktop import STKDesktop
    else:
        from agi.stk12.stkengine import STKEngine
    from agi.stk12.stkobjects import *
    from agi.stk12.stkobjects.astrogator import *
    from agi.stk12.utilities.colors import *
    from agi.stk12.stkengine.tkcontrols import GlobeControl, GfxAnalysisControl
except:
    print("Failed to import stk modules. Make sure you have installed the STK Python API wheel (agi.stk<..ver..>-py3-none-any.whl) from the STK Install bin directory")


if os.name == "nt":
    #获取根对象
    app = STKDesktop.StartApplication(visible=True, userControl=True)
    stkRoot = app.Root
else:
    app = STKEngine.StartApplication(noGraphics=False)
    stkRoot = app.NewObjectRoot()

if not stkRoot.AvailableFeatures.IsPropagatorTypeAvailable(AgEVePropagatorType.ePropagatorAstrogator):
    print("You do not have the required license.")
#######
#设置场景
Scenario_begin_time = "5 Mar 2024 04:00:00.000"#场景开始时间
Scenario_end_time = "9 Mar 2024 04:00:00.000"#场景结束时间
Scenario_epoch_time = Scenario_begin_time#场景历元默认设置为场景开始时间
Animation_begin_time = Scenario_begin_time#动画开始时间默认设置为场景开始时间
Animation_step_value = 10#动画步长
#新建一个场景
stkRoot.NewScenario('STK12')
#返回一个场景对象,以修改场景时间及其他属性
oScenario = stkRoot.CurrentScenario
oScenario.SetTimePeriod(Scenario_begin_time,Scenario_end_time)#设置场景时间间隔
oScenario.Epoch = Scenario_epoch_time#设置场景历元
#场景是场景,动画是动画,动画的开始时间也要单独设置
oScenario.Animation.StartTime = Animation_begin_time#设置为场景开始时间
oScenario.Animation.AnimStepValue = Animation_step_value#设置动画步长
stkRoot.Rewind()  # 将场景返回时间起点
#######
#读取TLE数据,构建卫星
# 地球引力常数(km^3/s^2)
GM = 398600.4418
#打开TLE数据文件
f = open('TLEfitting20200901_Simple.txt', 'r')
line = f.readlines()
for i in range(0, len(line), 2):
    #######
    #这段代码来自AI,读取TLE数据,得到卫星轨道六要素
    line1 = line[i].strip('\n')
    line2 = line[i + 1].strip('\n')
    satellite = EarthSatellite(line1, line2)
    # 获取卫星编号
    satellite_number = satellite.model.satnum
    if satellite_number<10:
        satellite_name = '0000'+str(satellite_number)
    elif 10<=satellite_number<100:
        satellite_name = '000'+str(satellite_number)
    elif 100<=satellite_number<1000:
        satellite_name = '00'+str(satellite_number)
    elif 1000<=satellite_number<10000:
        satellite_name = '0'+str(satellite_number)
    else:
        satellite_name = str(satellite_number)
    # 获取轨道参数
    model = satellite.model
    # 访问轨道参数
    inclination = math.degrees(model.inclo)  # 轨道倾角,单位°
    eccentricity = model.ecco#轨道偏心率
    right_ascension = math.degrees(model.nodeo)  # 升交点赤经,单位°
    argument_of_perigee = math.degrees(model.argpo)#近地点辐角
    mean_motion = model.no  # 平均飞行速度,弧度/分钟
    mean_anomaly = math.degrees(model.mo)#平近点角
    # 转换平均运动为每天圈数
    mean_motion_revs_per_day = mean_motion * 1440 / (2 * math.pi)
    # 计算半长轴(单位:km)
    mean_motion_rad_per_sec = mean_motion / 60  # 转换为弧度/秒
    semi_major_axis = (GM / (mean_motion_rad_per_sec ** 2)) ** (1 / 3)
    #######
    #新建一个卫星对象
    oSat_VI = oScenario.Children.New(AgESTKObjectType.eSatellite, satellite_name)
    oSat_VI.VO.Model.ScaleValue = 0.0
    oSat_VI.SetPropagatorType(AgEVePropagatorType.ePropagatorJ2Perturbation)  # 设置轨道传播器摄动类型J2
    oJ2Perturbation_VI = oSat_VI.Propagator  # 获取传播器对象(后面要设置属性)
    interval_VI = oJ2Perturbation_VI.EphemerisInterval  # 获取传播时间对象
    interval_VI.SetExplicitInterval(Scenario_begin_time, Scenario_end_time)  # 设置传播时间
    oJ2Perturbation_VI.Step = 60  # 设置计算步长
    #输入轨道六要素,预报轨道
    oOrb_VI = oJ2Perturbation_VI.InitialState.Representation.AssignClassical(3, semi_major_axis,
                                                                             eccentricity, inclination,
                                                                             argument_of_perigee, right_ascension,
                                                                             mean_anomaly)
    oJ2Perturbation_VI.Propagate()  # 传播计算
    oSat_VI.Graphics.SetAttributesType(AgEVeGfxAttributes.eAttributesBasic)
    SatVI_attributes = oSat_VI.Graphics.Attributes
    SatVI_attributes.Inherit = False
    SatVI_attributes.IsOrbitVisible = False
f.close()
#######
#为每颗卫星设置传感器,传感器信息文件可以自己设置
Sat_ele = oScenario.Children.GetElements(AgESTKObjectType.eSatellite)
sat_num = Sat_ele.Count
Sensor_Workbook = openpyxl.load_workbook('传感器信息.xlsx')
Sensor_Sheet = Sensor_Workbook.active
Sensor_data = []
Sensor_name = []
for Sensor_row in Sensor_Sheet.iter_rows(values_only=True):
    Sensor_data.append(Sensor_row)
del Sensor_data[0]
for i in range(len(Sensor_data)):
    Sensor_name.append(Sensor_data[i][0])
for i in range(sat_num):
    #获取卫星对象
    Sat_item = Sat_ele.Item(j)
    satPath = Sat_item.Path
    sat = stkRoot.GetObjectFromPath(satPath)
    sat_name = sat.InstanceName
    if sat_name in Sensor_name:
        x = Sensor_name.index(sat_name)
        if Sensor_data[x][1]=='O':
            satSen = oSat_VI.Children.New(20, sat_name+'O')
            satSen.SetPatternType(AgESnPattern.eSnRectangular)  # 设置传感器类型(矩形)
            satSen.Pattern.HorizontalHalfAngle = Sensor_data[x][3]
            satSen.Pattern.VerticalHalfAngle = Sensor_data[x][4]
            satSen.Pattern.AngularResolution = Sensor_data[x][2]
            satSen_2D = satSen.Graphics
            satSen_2D.Color = Colors.Blue
        elif Sensor_data[x][1]=='S':
            satSen = oSat_VI.Children.New(20, sat_name+'S')
            satSen.SetPatternType(AgESnPattern.eSnRectangular)  # 设置传感器类型(矩形)
            satSen.Pattern.HorizontalHalfAngle = Sensor_data[x][3]
            satSen.Pattern.VerticalHalfAngle = Sensor_data[x][4]
            satSen.Pattern.AngularResolution = Sensor_data[x][2]
            satSen_2D = satSen.Graphics
            satSen_2D.Color = Colors.Red

评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值