STK Starlink Instances

这里我们给几个社区提供的 Starlink 可视化平台:

• Starlink Coverage Tracker
  • Starlink的可视化集大成者，超神 👍
• Live Starlink Satellite
  • Starlink / OneWeb / GPS 的最直观体现，便于比较 🔥

基础设施构建与连接

效果展示

github repo

这个例子展示 基础星链构建、地面站构建、地面站与LEO卫星的直线连接

结果演示:

控制台输出:

代码解析

库引入

# 时间处理、操作系统相关
import os
import platform
import time

# STK的对象和枚举类型
from agi.stk12.stkobjects import (
    AgEClassicalLocation,
    AgEClassicalSizeShape,
    AgECvBounds,
    AgECvResolution,
    AgEFmCompute,
    AgEFmDefinitionType,
    AgEOrientationAscNode,
    AgESTKObjectType,
    AgEVePropagatorType,
)

# 轨道状态
from agi.stk12.stkutil import AgEOrbitStateType

启动STK引擎或桌面应用

startTime = time.time()

if platform.system() == "Linux":
    useStkEngine = True
else:
    useStkEngine = False

'''
根据操作系统的不同，代码选择启动STK引擎（在Linux上）或STK桌面应用（在Windows上）。useStkEngine变量用于控制使用哪种模式。
- 启动后，获取根对象stkRoot，这是所有其他STK对象的父级
'''

if useStkEngine:
    from agi.stk12.stkengine import STKEngine
    stk = STKEngine.StartApplication(noGraphics=True)
    stkRoot = stk.NewObjectRoot()
else:
    from agi.stk12.stkdesktop import STKDesktop
    stk = STKDesktop.StartApplication(visible=True, userControl=True)
    stkRoot = stk.Root

设置日期格式和创建场景

# 日期格式为UTC
stkRoot.UnitPreferences.SetCurrentUnit("DateFormat", "UTCG")
# 打开系统并新建Scenario
stkRoot.NewScenario("PythonEngineExample")
# 选定当前的Scenario
scenario = stkRoot.CurrentScenario
# 设定观测/运行对应的时间区间
scenario.SetTimePeriod("1 Aug 2020 16:00:00", "2 Aug 2020 16:00:00")

创建卫星并设置轨道参数

# 创建一个新的卫星对象
# scenario.Children.New(type, name): 在场景的子对象中添加新对象
# type: AgESTKObjectType.eSatellite
# name: MySatellite
satellite = scenario.Children.New(AgESTKObjectType.eSatellite, "MySatellite")

# 设置传播器类型
# SetPropagatorType(type)
# type: AgEVePropagatorType.ePropagatorTwoBody (两体传播器)
satellite.SetPropagatorType(AgEVePropagatorType.ePropagatorTwoBody)

# 获取传播器实例
propagator = satellite.Propagator

# 从传播器中获取初始轨道状态的表示 (location / velocity)
orbitState = propagator.InitialState.Representation
# 转换为经典轨道状态 (更适合进行传统的轨道分析和计算)
orbitStateClassical = orbitState.ConvertTo(AgEOrbitStateType.eOrbitStateClassical)

'''
这部分是轨道参数的设置，纯地理学名词，暂略
'''

# Set SMA and eccentricity
orbitStateClassical.SizeShapeType = AgEClassicalSizeShape.eSizeShapeSemimajorAxis
sizeShape = orbitStateClassical.SizeShape
sizeShape.Eccentricity = 0
sizeShape.SemiMajorAxis = 8000

# Set inclination and argument of perigee
orientation = orbitStateClassical.Orientation
orientation.Inclination = 25
orientation.ArgOfPerigee = 0

# Set RAAN and true anomaly
orientation.AscNodeType = AgEOrientationAscNode.eAscNodeRAAN
raan = orientation.AscNode
raan.Value = 0

orbitStateClassical.LocationType = AgEClassicalLocation.eLocationTrueAnomaly
trueAnomaly = orbitStateClassical.Location
trueAnomaly.Value = 0

# 导入轨道状态到传播器
# 传播器传播
# Assign orbit state and propagate satellite
orbitState.Assign(orbitStateClassical)
propagator.Propagate()

Note

1. 创建一个新的卫星对象
2. 初始化传播器
3. 获取传播器实例
4. 根据实例中的数据初始化轨道
5. 自定义参数表示轨道
6. 将轨道状态导入传播器
7. 传播器传播状态

创建设施并计算访问权限

# 在当前场景中创建一个新的设施对象，名称为 "MyFacility"
facility = scenario.Children.New(AgESTKObjectType.eFacility, "MyFacility")
# 设置设施的地理位置 (纬，经，高)
facility.Position.AssignGeodetic(28.62, -80.62, 0.03)

# 计算卫星与刚刚创建的设施之间的访问权限
# 返回值 access 包含了“访问计算”的相关信息
access = satellite.GetAccessToObject(facility)
# 根据卫星的轨道和设施的位置，确定卫星在给定时间段内是否能够“看到”该设施
access.ComputeAccess()
# 建立“地面站”与对应“LEO卫星”的链路连接

# Get access interval data and print results
# 提取并打印在指定时间范围内的所有有效访问间隔
stkRoot.UnitPreferences.SetCurrentUnit("Time", "Min")
accessDataProvider = access.DataProviders.GetDataPrvIntervalFromPath("Access Data")
elements = ["Start Time", "Stop Time", "Duration"]
accessResults = accessDataProvider.ExecElements(scenario.StartTime, scenario.StopTime, elements)

打印访问结果到控制台

# 获取 StartTime / StopTime / Duration
startTimes = accessResults.DataSets.GetDataSetByName("Start Time").GetValues()
stopTimes = accessResults.DataSets.GetDataSetByName("Stop Time").GetValues()
durations = accessResults.DataSets.GetDataSetByName("Duration").GetValues()

print("\nAccess Intervals")
for i in range(len(startTimes)):
    print("{a:<29s}  {b:<29s}  {c:<4.2f}".format(a=startTimes[i], b=stopTimes[i], c=durations[i]))

print("\nThe maximum access duration is {a:4.2f} minutes.".format(a=max(durations)))

对应结果

Access Intervals
Start Time                     Stop Time                      Duration (min)
1 Aug 2020 16:25:03.475805187  1 Aug 2020 16:51:07.039122204  26.06
1 Aug 2020 18:33:25.466968796  1 Aug 2020 18:57:49.382137265  24.40
1 Aug 2020 20:43:39.293864757  1 Aug 2020 21:02:22.893522448  18.73
2 Aug 2020 07:42:48.464807220  2 Aug 2020 07:58:38.079188917  15.83
2 Aug 2020 09:46:27.097228442  2 Aug 2020 10:09:58.224122333  23.52
2 Aug 2020 11:52:44.502495801  2 Aug 2020 12:18:34.702365047  25.84
2 Aug 2020 14:00:19.322304548  2 Aug 2020 14:26:38.819331701  26.32

The maximum access duration is 26.32 minutes.
--- Access computation: 6.141 sec    Total time: 64.791 sec ---

创建星座与链对象并计算访问权限

# 卸载初始卫星
satellite.Unload()
# 创建一个名为 “SatConstellation” 的星座对象
constellation = scenario.Children.New(AgESTKObjectType.eConstellation, "SatConstellation")

# ... (4 orbit, each with 8 satellites)
# 循环创建多个卫星并将其添加到星座中

# 创建链对象，将星座与设施添加到链中，并计算访问权限
# Insert the constellation of Satellites with specific parameters...
chain.ComputeAccess()

创建覆盖对象并计算覆盖率

# 创建覆盖对象
coverageDefinition = scenario.Children.New(AgESTKObjectType.eCoverageDefinition, "CoverageDefinition")
grid.BoundsType = AgECvBounds.eBoundsCustomRegions 

# Add US shapefile to bounds
# 边界是USA (自输入对应的文件)
bounds = coverageDefinition.Grid.Bounds

if platform.system() == "Linux":
    install_path = os.getenv("STK_INSTALL_DIR")
else:
    import winreg

    registry = winreg.ConnectRegistry(None, winreg.HKEY_LOCAL_MACHINE)
    key = winreg.OpenKey(registry, r"Software\AGI\STK\12.0")
    install_path = winreg.QueryValueEx(key, "InstallHome")

bounds.RegionFiles.Add(
    os.path.join(
        install_path[0],
        r"Data/Shapefiles/Countries/United_States_of_America\United_States_of_America.shp",
    )
)

# set resolution...
coverageDefinition.ComputeAccesses()

可视化效果 (图中紫色区域, 即美国国界)

输出覆盖率计算结果并关闭STK

minAccessDuration = fomResults.DataSets.GetDataSetByName("Minimum").GetValues()[0]
maxAccessDuration = fomResults.DataSets.GetDataSetByName("Maximum").GetValues()[0]
avgAccessDuration = fomResults.DataSets.GetDataSetByName("Average").GetValues()[0]

print("\nThe minimum coverage duration is {a:4.2f} min.".format(a=minAccessDuration))
print("The maximum coverage duration is {a:4.2f} min.".format(a=maxAccessDuration))
print("The average coverage duration is {a:4.2f} min.".format(a=avgAccessDuration))

stkRoot.CloseScenario()
stk.ShutDown()
print("\nClosed STK successfully.")

对应结果

The minimum coverage duration is 20.91 min.
The maximum coverage duration is 24.43 min.
The average coverage duration is 23.65 min.
--- Coverage computation: 17.676 sec    Total time: 96.185 sec ---

What's More

1) 如果我们在现在结果的基础上，加上根据TLE生成的Starlink，会如何？

结果清晰地说明了这两者会 同时显示在当前Scenario上，互不影响

2) 之前的配置是 4 orbit. 8 sat/orbit. 现在改成 10 orbit. 10 sat/orbit.

我们清晰地看出， 轨道数量、单轨卫星数、地面站连接覆盖的卫星数 都发生了相应的变化

Works Well 🎉 👍

TL; DR

• 创建卫星并设置轨道参数
• 创建地面站，建立 GS-LEO 链路
• 创建星座与链, 计算访问权限
• 创建覆盖对象, 计算指定区域覆盖率

代码生成的 和 文件导入的 互不影响

自定义 Starlink 星系

建立一个16x16的 starlink 星座拓扑

• 16 orbit shell
• each shell has 16 satellites

import time
# Create a Progress Bar on Console
from tqdm import tqdm
# STK Related
from comtypes.gen import STKObjects, STKUtil, AgSTKVgtLib
# Object Related
from comtypes.client import CreateObject, GetActiveObject, GetEvents, CoGetObject, ShowEvents
# Comtypes Related
from ctypes import *
from comtypes import GUID
from comtypes import helpstring
from comtypes import COMMETHOD
from comtypes import dispid
from ctypes.wintypes import VARIANT_BOOL
from ctypes import HRESULT
from comtypes import BSTR
from comtypes.automation import VARIANT
from comtypes.automation import _midlSAFEARRAY
from comtypes import CoClass
from comtypes import IUnknown
from comtypes.automation import IDispatch
from comtypes import DISPMETHOD, DISPPROPERTY, helpstring

############################################################################
# Scenario Setup
############################################################################

startTime = time.time()
useSTKEngine = False
readScenario = False

if useSTKEngine:
    # For Linux
    # (1) launch STK Engine
    print('Launching to STK Engine (STK12.2) ...')
    stkApp = GetActiveObject('STK12.Application')
    # (2) Disable graphics (on linux)
    stkApp.NoGraphics = True
    # (3) Ready for a new scenario
    print('Making a new stkRoot ...')
    stkRoot = CreateObject('AgStkObjects12.AgStkObjectRoot')
else:
    # For Windows
    # (1) launch STK
    print('Launching to STK (STK12.2) ...')
    if not readScenario:
        # creating a new scenario
        uiApp = CreateObject('STK12.Application')
    else:
        # from an existing scenario
        uiApp = GetActiveObject('STK12.Application')
    uiApp.Visible = True # GUI
    uiApp.UserControl = True # Mouse and GUI Console
    # (2) Ready for a new scenario
    print('Making a new stkRoot ...')
    stkRoot = uiApp.Personality2

'''
When creating a new scenario, the following code is needed
# For Linux (no GUI, create a new one):
stkRoot = CreateObject('AgStkObjects12.AgStkObjectRoot')
# For Windows (has GUI, complete instance):
stkRoot = stkApp.Personality2
'''

############################################################################
# Scenario Init
############################################################################

stkRoot.UnitPreferences.SetCurrentUnit('DateFormat', 'UTCG')
# Create a new scenario
print("Creating a new scenario ...")
if not readScenario:
    stkRoot.NewScenario('StarLink')

# eScenario Obj
scenario = stkRoot.CurrentScenario # python obj
print(scenario)
# IAgScenario Obj
scenario2 = scenario.QueryInterface(STKObjects.IAgScenario) # COM obj
print(scenario)

timeDuration = time.time() - startTime
splitTime = time.time()

print("--- Scenario creation: {a:4.3f} "
    "sec\t\tTotal time: {b:4.3f} sec ---".format(a=timeDuration, b=timeDuration))

'''
eScenario and IAgScenario:

- eScenario is a Python obj (object)
- IAgScenario is a COM obj (interface)

!!! You can only use IAgScenario to access STK Objects !!!
'''

############################################################################
# Constellation Setup
############################################################################

def CreateSatellite(
    numOrbitPlane=16, 
    numSatPerPlane=16, 
    hight=550, 
    inc=53,
    name="sat"
):
    # Create constellation object
    # eConstellation Obj
    constellation = scenario.Children.New(STKObjects.eConstellation, name)
    # IAgConstellation Obj
    constellation2 = constellation.QueryInterface(STKObjects.IAgConstellation)

    # Insert sat on each orbit plane
    for orbitIdx in range(numOrbitPlane):
        for satIdx in range(numSatPerPlane):
            # (1) Init a new satellite
            # eSatellite Obj
            satellite = scenario.Children.New(STKObjects.eSatellite, 
                                            f"{name}{orbitIdx}-{satIdx}")
            # IAgSatellite Obj
            satellite2 = satellite.QueryInterface(STKObjects.IAgSatellite)

            # (2) Two-body Propagator: IAgSatellite -> eSatellite
            satellite2.SetPropagatorType(STKObjects.ePropagatorTwoBody)
            twoBodyPropagator = satellite2.Propagator.QueryInterface(STKObjects.IAgVePropagatorTwoBody)

            # (3) Init state
            keplarian = twoBodyPropagator.InitialState.Representation.ConvertTo(
                STKUtil.eOrbitStateClassical).QueryInterface(STKObjects.IAgOrbitStateClassical)
            # a:
            keplarian.SizeShapeType = STKObjects.eSizeShapeSemimajorAxis
            # h:
            keplarian.SizeShape.QueryInterface(
                STKObjects.IAgClassicalSizeShapeSemimajorAxis).SemiMajorAxis = hight + 6371
            # e:
            keplarian.SizeShape.QueryInterface(
                STKObjects.IAgClassicalSizeShapeSemimajorAxis).Eccentricity = 0
            # inclination:
            keplarian.Orientation.Inclination = int(inc)
            # Arg of Perigee (近地点):
            keplarian.Orientation.ArgOfPerigee = 0
            # 升交点:
            keplarian.Orientation.AscNodeType = STKObjects.eAscNodeRAAN
            # 使卫星均匀分布在不同轨道面上
            RAAN = 360 / numOrbitPlane * orbitIdx
            keplarian.Orientation.AscNode.QueryInterface(
                STKObjects.IAgOrientationAscNodeRAAN).Value = RAAN
            # 真近点角:
            keplarian.LocationType = STKObjects.eLocationTrueAnomaly
            # 使卫星在同一轨道面上均匀分布
            trueAnomaly = 360 / numSatPerPlane * satIdx
            keplarian.Location.QueryInterface(STKObjects.IAgClassicalLocationTrueAnomaly).Value = trueAnomaly

            # (4) Propagate
            # All args -> keplarian -> twoBodyPropagator
            twoBodyPropagator.InitialState.Representation.Assign(keplarian)
            # Now do it
            twoBodyPropagator.Propagate()

            # (5) Add to constellation
            constellation2.Objects.AddObject(satellite)

############################################################################
# Sender and Receiver for Each Satellite
############################################################################

# Add Sender and Receiver for each Satellite
def AddTransmitterReceiver(sat_list):
    for each in sat_list:
        name = each.InstanceName
        # new transmitter and receiver
        transmitter = each.Children.New(STKObjects.eTransmitter, "Transmitter_" + name)
        reciver = each.Children.New(STKObjects.eReceiver, "Reciver_" + name)

# Sender Args
def SetTransmitterParameter(
    transmitter,
    frequency=12, 
    EIRP=20, 
    DataRate=14
):
    # 建立发射机的映射，以便对其进行设置
    transmitter2 = transmitter.QueryInterface(STKObjects.IAgTransmitter)
    # 选择 简单发射机 模型
    transmitter2.SetModel('Simple Transmitter Model')
    # Get Sending Model
    txModel = transmitter2.Model # py obj
    txModel = txModel.QueryInterface(STKObjects.IAgTransmitterModelSimple) # COM obj
    # Args
    txModel.Frequency = frequency
    txModel.EIRP = EIRP
    txModel.DataRate = DataRate

# Receiver Args
def SetReceiverParameter(receiver, GT=20, frequency=12):
    # 建立接收机的映射，以便对其进行设置
    receiver2 = receiver.QueryInterface(STKObjects.IAgReceiver)
    # 选择 简单接收机 模型
    receiver2.SetModel('Simple Receiver Model')
    # Get Sending Model
    recModel = receiver2.Model # py obj
    recModel = recModel.QueryInterface(STKObjects.IAgReceiverModelSimple) # COM obj
    # Args
    recModel.AutoTrackFrequency = False
    recModel.Frequency = frequency
    recModel.GOverT = GT
    return receiver2

# Get Receiver Instances
def GetSatReceiver(sat, GT=20, frequency=12):
    # py obj
    receiver = sat.Children.GetElements(STKObjects.eReceiver)[0]
    # COM obj
    receiver2 = SetReceiverParameter(receiver=receiver, GT=GT, frequency=frequency)

    return receiver2

############################################################################
# Compute Access (link-level)
############################################################################

def LinkComputeAccess(access):
    # Init for Access
    access.ComputeAccess()
    # Get DataProvider
    accessDP = access.DataProviders.Item('Link Information')
    accessDP2 = accessDP.QueryInterface(STKObjects.IAgDataPrvTimeVar)

    Elements = ["Time", 'Link Name', 'EIRP', 'Prop Loss', 'Rcvr Gain', "Xmtr Gain", "Eb/No", "BER"]
    results = accessDP2.ExecElements(scenario2.StartTime, scenario2.StopTime, 3600, Elements)

    Times = results.DataSets.GetDataSetByName('Time').GetValues()
    EbNo = results.DataSets.GetDataSetByName('Eb/No').GetValues()
    BER = results.DataSets.GetDataSetByName('BER').GetValues()
    Link_Name = results.DataSets.GetDataSetByName('Link Name').GetValues()
    Prop_Loss = results.DataSets.GetDataSetByName('Prop Loss').GetValues()
    Xmtr_Gain = results.DataSets.GetDataSetByName('Xmtr Gain').GetValues()
    EIRP = results.DataSets.GetDataSetByName('EIRP').GetValues()

    return Times, Link_Name, BER, EbNo, Prop_Loss, Xmtr_Gain, EIRP

# Get all Access
accessList = scenario2.GetExistingAccesses()
for accessPath in tqdm(accessList):
    # accessPath: ('Satellite/Sat0_0/Transmitter/Transmitter_Sat0_0','Satellite/Sat0_1/Receiver/Reciver_Sat0_1',True)
    # accessPath[0]: sender address
    # accessPath[1]: receiver address
    transmitterName = accessPath[0].split('/')[-1]
    reciverName = accessPath[1].split('/')[-1]
    access = scenario2.GetAccessBetweenObjectsByPath(accessPath[0], accessPath[1])
    link = access
    transmitterName = accessPath[0].split('/')[-1]
    reciverName = accessPath[1].split('/')[-1]
    access = scenario2.GetAccessBetweenObjectsByPath(accessPath[0], accessPath[1])

############################################################################
# Color the Satellite
############################################################################

def Change_Sat_color(sat_list):
    print('Changing Color of Satellite ...')

    for each_sat in tqdm(sat_list):
        now_sat_name = each_sat.InstanceName
        now_plane_num = int(now_sat_name.split('_')[0][3:])
        now_sat_num = int(now_sat_name.split('_')[1])

        satellite = each_sat.QueryInterface(STKObjects.IAgSatellite)
        graphics = satellite.Graphics
        graphics.SetAttributesType(1)  # eAttributesBasic
        attributes = graphics.Attributes
        attributes_color = attributes.QueryInterface(STKObjects.IAgVeGfxAttributesBasic)
        attributes_color.Inherit = False

        color_sheet = [16436871, 2330219, 42495, 9234160, 65535, 255, 16776960]
        if now_sat_name[2] == 'A':
            color = 255
        elif now_sat_name[2] == 'B':
            color = 42495
        elif now_sat_name[2] == 'C':
            color = 16436871
        attributes_color.Color = color
        # orbit attribute interface
        orbit = attributes.QueryInterface(STKObjects.IAgVeGfxAttributesOrbit)
        orbit.IsOrbitVisible = False

############################################################################
# Main
############################################################################

if not readScenario:
    CreateSatellite(numOrbitPlane=16, numSatPerPlane=16, hight=550, inc=53)
    sat_list = stkRoot.CurrentScenario.Children.GetElements(STKObjects.eSatellite)
    AddTransmitterReceiver(sat_list)

    # Change DateFormat dimension to epoch seconds to make the data easier to handle in
    # Python
    stkRoot.UnitPreferences.Item('DateFormat').SetCurrentUnit('EpSec')
    # Get the current scenario
    scenario = stkRoot.CurrentScenario
    # Set up the access object
    access = constellation2.GetAccessToObject()
    access.ComputeAccess()
    # Get the Access AER Data Provider
    accessDP = access.DataProviders.Item('Access Data').Exec(scenario.StartTime, scenario.StopTime)

    accessStartTimes = accessDP.DataSets.GetDataSetByName('Start Time').GetValues
    accessStopTimes = accessDP.DataSets.GetDataSetByName('Stop Time').GetValues

效果演示 (正视 + 从北极俯视)