CN102841966A - Vpp-STK卫星仿真开发与运行平台系统 - Google Patents

Abstract

发明名称：Vpp-STK卫星仿真运行与平发平台系统。所属的技术领域：电子信息-计算机及网络产品-空间信息获取及综合应用集成系统。所要解决的技术问题：本发明主要以国产技术研发的专利产品，替代美国长期对中国技术封锁限制的美国Analytical Graphics公司开发的STK(Satellite Tool Kit)卫星工具软件、打破美国对中国禁售、禁运的封锁。解决该问题的技术方案的要点：本发明基于国产VBF组件化框架设计、底层基于C/C++研制开发，做为中国航天领域进行卫星仿真科研论证支撑平台。主要用途：实现卫星、运载火箭、导弹、地面站的飞行器卫星轨道分析与数据计算并可在二维地图和三维显示环境显示。

Description

Vpp-STK卫星仿真开发与运行平台系统

一、技术领域

所属的技术领域：电子信息-计算机及网络产品-空间信息获取及综合应用集成系统。 

本发明专利是一种用于卫星仿真开发与运行的工具、是国产的卫星仿真开发与运行平台。 

本发明专利用于打破美国长期以来对中国禁售禁运的卫星仿真技术(美国Analytical Graphics公司开发的STK(Satellite Tool Kit)卫星工具软件)，解决国产化、替代美国商业仿真软件产品的目的。主要应用技术领域： 

①卫星、运载火箭、导弹、飞机、地面目标等生成位置和姿态数据、可见性及遥控器覆盖分析提供分析引擎，用于任务数据计算分析并可在二维地图和三维显示环境显示。 

②卫星、运载火箭、导弹的轨道预报算法、姿态定义、坐标类型和坐标系统、遥感器类型，以及卫星、城市、地面站和恒星数据库。 

③卫星、运载火箭、导弹、地面站之间的通信分析、雷达分析、覆盖分析、轨道机动、精确定轨、实时操作等高级约束条件定义对于特定的分析。 

二、背景技术

目前国内使用的是一款美国的卫星仿真开发工具，而这款工具自2004年以来一直对华技术封锁、禁售、禁运。 

本发明专利是就是一种国产平台的、用于卫星仿真开发与运行的工具。打破美国对华的技术封锁与垄断。目前国内尚无卫星仿真商业软件、工具、平台。具体的技术背景如下所述： 

一是美国FBI对华禁售、禁运：目前美国Analytical Graphics公司开发的STK(Satellite Tool Kit)卫星工具软件(以下简称美国STK)，从2004年的STK v4.0版开始对中国禁售、禁运、技术封锁至今(目前发展到Ver9.2版)。封锁对象包括美国STK分为基本版(对应应用技术领域1)、专业版(对应应用技术领域2)、高级版(对应应用技术领域3)。 

虽经费与实力强的航天领域用户则从第三方国家采购到，正版软件(对应应用技术领域2)，但是无法得到正版软件的应该有的技术支持与升级服务。如美国STK不提供中文显示支持、也不可能提供中国售后服务。 

三是受制于美国对中国的技术封锁：航天研究相关部门的开发，受限于美国对中国的STK技术封锁限制，特别是高级版无法从第三国采购到(对应应用技术领域3)。如果进行深入的研究、修改完善其关键点内容，则无法在底层源代码的深度上满足我国国防军工研究使用，且不支持中国54/80/CGCS2000三大体系的坐标系投影标准。 

四是国产化打破美国封锁：为打破美国对中国的STK卫星仿真技术的封锁，自2005年开始研究美国的STK的技术结构与功能组成，根据中国航天工业科研的实际，设计研制了一套底层基于C/C++语言开发的纯国产的Vpp-STK卫星仿真运行与平发平台工具，既可以满足中国的航天用户使用，又能满足宇航专业的大学和研究所使用，同时还能在国产linux操作系统和国产嵌入式操作系统VxWorks上运行，彻底解决了跨操作系统平台、中文显示、支持中国标准的坐标投影换算、卫星轨道与姿态的高精度解算等一系列难题。 

五是回避了美国过旧的技术垃圾：美国NASA(National Aeronautics and Space Administration)美国国家航空及太空总署，将技术已经过时的WorldWind.exe源代码开放到网上，供天文爱好者研究讨论；国内一些研究所和遥感影像数据代理公司借此机会下载，并基于NASA这套开源代码、做了一些学习研究；但真正高精尖技术是不可能放在网上供外人、特别是中国人研究掌握和消化的。 

另外，美国FBI默许美国学术界的6名研究人员将OSG(Open Scene Graph)源代码放在网上，供爱好者研究、开发，根据美国对华高技术封锁一惯示例，有价值的技术特别是源代码是根本不可能放在网上，放在网上发布的OSG和NASA开源代码，都不包含美国STK卫星仿真软件的关键技术与工具。 

即使使用了WorldWind和OSG，因是纯学术的、探讨研究性质的，基于BBS论坛发布博客、微博进行交流经验与体会的，是坊间的学术性质的、不成熟的半成品，航天用户无法在非成品上进行工程实用，也无法得到所需要的技术支持与售后服务，因为开源的代码是不需要负责任、也不需要负责售后的，更不是商业软件工具、不可能包含美国的核心技术精华内容。 

三、发明内容

本发明专利是一种国产平台的、用于卫星仿真开发与运行的工具，属于商业化的卫星仿真平台系统。具体内容如下： 

(1)本发明专利是以国产技术研发的“Vpp-STK卫星仿真运行与平发平台系统”专利产品，替代美国长期对中国技术封锁与技术限制控制的、航天卫星仿真开发的STK(Satellite Tool Kit卫星工具包)，打破美国对中国禁售、禁运的封锁。 

(2)Vpp-STK卫星仿真运行与平发平台系统，底层基于C++、使用VVP-VBF组件化开发框架开发，解决了中文汉字显示、底层平台不开放、跨操作系统平台移植、对中国禁售禁运、技术开发服务不支持等一系列问题。 

(3)在国产Vpp-STK卫星仿真运行与平发平台系统上，实现飞行器卫星轨道分析与三维空间任务显示、卫星探测覆盖分析、预警探测轨道分析、接近碰撞分析、拦截器飞行分析、来袭目标弹道导弹 分析、链路分析、通信分析、电磁分析、雷达分析、推演分析。 

(4)在国产Vpp-STK卫星仿真运行与平发平台系统上，实现卫星、地面站和日、月、地、行星、恒星星历库加载显示。支持查询卫星相关参数，包括轨道类型、轨道参数和卫星属性、具体位置信息，卫星编号、时间点、各时间点位置坐标、各时间点速度矢量等。 

(5)在国产Vpp-STK卫星仿真运行与平发平台系统上，实现地面站测控分析，支持坐标系换算(地心惯性坐标系、地心固定坐标系、地心轨道坐标系、黄道坐标系、卫星发射坐标系，及其空间配准)，默认为地心惯性坐标系，支持时间系统可选(太阳时、恒星时、国际原子时、动力学时，默认为格林尼治恒时)(目前提供世界时)。 

四、附图说明

参照以下说明、所附权利要求以及附图，将更好地理解本发明专利的功能和方法，这些设计构思及其技术特征，说明书附图中给出本发明专利的各功能模块之间的计算机程序的主要流程图，技术方案、整体设计、工作过程和步骤的技术特征、优点方面，在发明创新系统的附图中： 

图1是本发明专利的技术，在航天城使用本发明专利系统软件进行多颗卫星轨道运行、对地观测、覆盖分析仿真计算界面。该发明专利已实现国外卫星仿真软件STK的空间飞行任务与空间三维态势可视化环境显示效果；星空背景是行星的星历库文件加载后的效果。 

图2是本发明专利的技术，在航天704研究所使用该专利系统、实现近地空间的地球同步轨道通信卫星仿真、航天空间飞行器飞行仿真、空间飞行器地面测控站仿真、以及空间飞行器与地面站之间通信关系实时仿真分析的实际使用界面。 

图3是本发明专利的技术、基于国产Vpp-STK卫星仿真运行与平发平台系统的、在二维GIS地理信息系统上的全世界地图上的星下点轨迹生成图；达到国际卫星仿真的先进水平。 

图4是本发明专利的国产Vpp-STK卫星仿真运行与平发平台系统组成与模块组件、层次结构组成框图，是实现图1、图2和图3中功能的支撑。是对图5的技术架构实现的功能说明的详细描述。采用分层的体系结构，分为运行环境抽象层、核心层、界面层、COM组件封装层、体系结构适配层等五个层次，采用这样的体系结构带来的优势是： 

■跨平台的移植能力，可支持Windows、Linux、VxWorks等操作系统； 

■可裁剪能力，适用于多种硬件平台，特别是可应用在资源有限的嵌入式平台上； 

■支持多种应用模式，可以满足多种体系架构的需要； 

■灵活的二次开发模式，可在多个层次上进行二次开发和扩充平台功能； 

■很好的解决了多领域合作的问题。使得我们能把更多的精力集中于战场环境仿真专业领域的研究，而各种其他专业(如电子对抗)的模型由各专业的人员建立。 

图5是本发明专利的技术实现原理图，是对图4中的本发明专利的国产Vpp-STK卫星仿真运行与 平发平台系统技术原理与组成进行说明；是对图6的VBF组件化开发模式框架的专业描述。 

图6是本发明专利的所采用的VBF组件化开发模式框架的组成结构、系统结构、功能层次的说明，实现的依据、进行系统技术原理说明。 

图7是本发明专利的卫星轨道生成模型的人机交互界面，是对图4中卫星轨道分析功能模块的实际实现运行效果。 

图8是本发明专利的系统的星下点轨迹在三维空间态势上的飞行与显示，是对应图3中二维GIS的星下点轨迹生成图的三维空间显示。 

图9和图10是本发明专利的导弹弹道预测分析、拦截器飞行工具分析、接近碰撞分析、来袭目标(弹道导弹)分析、推演分析功能模块实现，是对应图4中系统功能描述的具体实现。 

图11和图12是本发明专利的链路分析、通信分析、电磁态势分析、雷达分析功能实现，以及在三维空间的仿真计算实时可视化。 

图13是本发明专利的数字勘场与全过程仿真应用，应用于神州八号、九号与天宫一号的近地空间飞行器测控仿真与全过程可视化的技术效果。 

图14是本发明专利的卫星轨道的生成算法组件的算法可视化，构造卫星轨道的经典算法，做为例程源代码供航天专业用户在此基础上进行算法模型的修正。 

图15是本发明专利的卫星轨道生成算法中的偏近点角和真近点角F的关系 

图16是本发明专利的卫星轨道生成算法中的牛顿法进行迭代求解的求解方法 

图17是本发明专利的卫星轨道生成算法中的生成一条二体轨道曲线流程 

图18是本发明专利的低轨卫星二体轨道惯性坐标系和地理坐标系下的生成轨道区别 

图19是本发明专利的卫星轨道生成的低轨卫星轨道的星点下点轨迹二维显示 

图20是本发明专利的卫星轨道生成的高轨卫星轨道(惯性坐标系和地理坐标系) 

图21是本发明专利的卫星轨道生成的高轨卫星轨道星下点轨迹 

图22是本发明专利的卫星轨道的卫星的地心赤道惯性坐标系下位置坐标在地理坐标系下的表示 

图23是本发明专利的卫星轨道的近地轨道地球形状J₂摄动 

图24是本发明专利的卫星轨道的近地轨道地球形状J₂摄动一个月以后的轨道。 

图25是本发明专利的无摄二体轨道动力学数值积分方法求解初始位置速度问题流程图(动力学积分模块)。

图26是本发明专利的解决摄动问题流程图。 

五、具体实施方式

本发明专利具体实现方式：国产Vpp-STK卫星仿真运行与平发平台系统 

(1)底层基于C/C++开发：是中国研制生产的、纯国产化的、基于国产化 

组件化开发框架结构的、底层用C/C++开发的、集成了一套卫星仿真中间件组件的国产Vpp-STK卫星仿真开发与运行平台。 

(2)先进的组件化框架结构：国产 

组件化开发框架系统内核，采用主控/组件与插件相结合的开放式总体架构，实现组件化开发模式，通过主控+注册组件+服务消息调用机制组件框架、与接口组件之间的消息发送、及接口查询，都必须通过主控完成组件注册、接口查询、命令发送、消息传递、事件监听、事件处理等，可实现多人多组协同组件化开发模式，大大提高了并行开发效率。这种软件结构可伸缩的技术架构，软件功能可随时扩展、大小结构可伸缩，便于卫星仿真不同类型的不同权限用户开发集成部署。 

(3)基于组件化开发向导模式：提供卫星仿真开发向导、组件式开发示例，实现卫星轨道分析与三维空间任务显示、覆盖分析、导弹弹道预测分析、拦截器飞行工具分析、接近碰撞分析、来袭目标(弹道导弹)分析、链路分析、通信分析、电磁态势分析、雷达分析、推演分析功能模块组件，通过注册组件形式、集成至Vpp-STK卫星仿真开发平台框架上，实现打破了美国对中国的卫星仿真技术封锁、航天技术限制。国产Vpp-STK卫星仿真运行与平发平台系统提供的功能，覆盖了国外(美国STK)的高级版主要功能，而这正是中国国防军工研究的关键领域。 

(4)跨操作系统平台：基于 

组件化开发框架结构的、基于主控+注册组件架构、底层基于C/C++开发的一套卫星仿真开发与运行平台，可移植到QT+liunx操作系统环境和QT+VxWorks嵌入式实时操作系统环境上运行，可接国产数据库、可移植到国产操作系统上，真正实现全国产化国防军工使用。 

支持目前主流开发工具MFC、BCG、C#、QT4.x和VC++6.0、VS2005/2008/2010、C/C++开发语言。支持windows XP/windows 7操作系统运行，嵌入式操作系统VxWorks5.5/6.6、Linux、国产红旗Linux操作系统的跨平台编译。 

(5)开放式集成框架模式：采用VVP-VBF组件化框架平台是能够与二维MGIS和三维GIS无缝集成，可集成国产VVP-3D三维视景仿真引擎v1.0、Deep Eye三维战场空天态势显示平台v3.0，可集成航天研究人员基于VC++6.0开发的成果。支持海量数据的快速显示：支持海量矢量数字地图数据、像素数字地图数据、数字高程模型数据、卫星遥感正射影像图数据加载显示，以及非系列比例尺地图(矢量形式或栅格形式)的检索与显示。 

(6)满足卫星仿真研究的二次开发：集成日、月、地、行星、恒星星历库，对太阳、地球、月亮三者间的相对运动关系进行分析，在地球上形成的天光地影及白天、黑夜等现象的时空关系，可在三维数字地球上进行日月地关系可视化显示。 

支持对地观测(扫描地球)、恒星观测、月球观测、冷空间观测等空间任务的多视角观测。支持地球观测实时数据和恒星实时数据(包括赤经、赤纬、亮度等)显示，能够配合数据驱动模块，产生生动逼真的三维空间探测任务显示。 

本发明专利的关键技术实现步骤： 

1、基于 

组件化开发框架结构 

(1)熟悉组件化开发框架VBF接口。VBF组件框架采用多层次主控/组件结构，由主控(只有一个)、二级主控(可能多个)及若干组件组成，控主要完成组件注册、消息分发、组件通信、接口查询等的功能。组件完成具体的业务，可以接收并处理主控传来的各种命令、消息，并提供功能接口并外部调用。二级主控既是主控的一个组件，又是其管理的组件的主控。见图6 

(2)掌握组件化开发模式：系统内核采用主控/组件与插件相结合的开放式总体架构。每个组件完成一个完整的业务，通过添加组件进行业务的扩展。每个插件点完成一个业务环节(一类算法)，通过增加插件进行算法的扩展。使用组件化开发，是采用主控+注册组件+动态库形式的、软件结构可伸缩的技术架构，可实现多人同时分模块开发、采用开发向导、生成标准的组件，可动态加载运行，能够大大减少重复开发工作量，提高代码的复用率、模块组件的复用率。 

(3)使用组件化开发向导，根据设计的各个卫星仿真业务流程，通过使用组件化开发向导，采用VC++6.0开发，实现卫星仿真各个专业功能组件，使用脚本化和图形化方式操作调用，完成卫星仿真概念设计与方预案任务验证。 

(4)快速集成卫星仿真平台：通过主控+注册组件的机制，集成使用组件化开发向导工具、自动代码生成的卫星仿真模型与分析工具组件、卫星仿真中间件组件，形成Vpp-STK卫星仿真开发与运行平台的一套卫星仿真开发包。见图4和图5 

2、采用先进的软件结构设计实现 

卫星仿真开发与运行平台采用分层的体系结构，分为运行环境抽象层、核心层、界面层、COM组件封装层、体系结构适配层等五个层次：(见说明书附图4、图5) 

(1)运行环境抽象层主要负责封装操作系统及硬件的细节，提高可移植性。同时，使得上层应用的开发可以在普通微机上进行，开发成本也因此大大降低。 

(2)核心层是与战场环境仿真专业相关的部分。该部分使用内建的对象与接口的注册、发现和执行机制，在数据模型支持、战场实体建模等多个方面都采用插件模式。 

(3)界面封装层是对核心层人机交互的封装，该层实现各种窗体、控件及交互方式。界面层是多语言的，支持MFC、BCG、C#、QT4.x和VC++6.0、VS2005/2008/2010并可实现语言之间的热切换。 

(4)COM组件封装层目的是实现对多种开发语言和平台的支持，支持MFC、BCG、C#、QT4.x和VC++6.0、VS2005/2008/2010，并充分利用COM+的各种优势，为组件进一步封装提供基础。 

(5)体系结构适配层使得VBF能够在各种体系结构下使用。如，采用传统体系结构的应用可以基于三个层次开发，一是基于核心层，采用C++语言开发，自己建立交互界面，使用内部提供的对象与接口的注册、发现与执行机制；二是基于界面封装层使用C++开发；三是基于COM组件封装层开发，基于该层可以使用支持COM的所有开发语言进行开发。 

采用这样的体系结构带来的优势是： 

■跨平台的移植能力，可支持Windows、Linux、VxWorks等操作系统； 

■可裁剪能力，适用于多种硬件平台，特别是可应用在资源有限的嵌入式平台上； 

■支持多种应用模式，可以满足多种体系架构的需要； 

■灵活的二次开发模式，可在多个层次上进行二次开发和扩充平台功能； 

很好的解决了多领域合作的问题。使得我们能把更多的精力集中于战场环境仿真专业领域的研究，而各种其他专业(如电子对抗)的模型由各专业的人员建立。 

3、以卫星轨道分析与三维空间任务显示功能模块组件为例，说明关键技术设计 

(1)功能描述 

能够提供卫星数据库(在轨、不工作卫星)轨道参数(建立一个卫星数据库用于记录入库，用于新建一个卫星、通过初始设置条件设置生成一个轨道)； 

能够对卫星等各类航天器等空间目标的轨道进行查询、编辑(增删改)； 

能够直接导入真实的卫星轨道数据模型文件。 

(2)性能要求 

●卫星数据库容量：支持20000条卫星等空间目标数据； 

●支持空间任务的多视角观测：能够配合其他模块数据驱动，产生生动、逼真的三维空间任务显示。 

(3)工作流程 

①卫星轨道生成及入库过程 

●选择轨道类型(太阳同步、对地静止、椭圆轨道、圆轨道) 

●选择轨道预报方法(二体，J2&J4，MSGP4，HPOP，LOP，外部文件) 

●选择卫星属性(未指定、敌、我、友、不明) 

●输入卫星轨道参数(半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角、过近地点时刻、高度、近地点、远地点、半径、偏心率) 

●设置卫星轨道颜色、线型 

●卫星轨道入库 

②卫星轨道查询过程 

●输入：待查询可选项：卫星编号、卫星轨道类型、卫星属性 

●输出：查询的卫星相关参数，包括轨道类型、轨道参数和卫星属性等 

③卫星轨道数据导入 

●卫星设置界面增加卫星轨道数据导入界面，通过导入的卫星轨道数据推出对应卫星轨道参数，并将卫星入库。 

●读轨道数据文件，此功能未完成，轨道数据文件格式不是轨道六根数，是具体位置信息，卫星编号、时间点、各时间点位置坐标、各时间点速度矢量 

(4)输入输出 

输入：轨道预报方法、轨道类型、卫星属性、卫星轨道参数、卫星轨道颜色、线型、设置卫星开始运行时刻、设置卫星初始化位置参数、卫星结束运行时刻 

输出：卫星轨道及各时间点卫星位置及速度矢量，二维和三维显示环境下显示卫星飞行效果。 

(5)约束条件 

●坐标系可选(地心惯性坐标系、地心固定坐标系、地心轨道坐标系)，默认为地心惯性坐标系 

●时间系统可选(太阳时、恒星时、国际原子时、动力学时)，默认为格林尼治恒星时(目前提供世界时) 

●卫星RCS留有可输入参数对话框，默认为0.1平方米 

●卫星姿态为标准姿态 

●卫星轨道数据文件格式 

●数据精确度小数点后6位 

(6)卫星轨道的生成算法 

卫星无摄二体轨道以万有引力公式为基础，假设地球和卫星都是理想质点的基础上得到的理想轨道，其表示方法有许多种，但是都以轨道六个数为基础，给出轨道六根数，就能得到卫星轨道。 

a    轨道半长轴，单位：km 

e    轨道偏心率 

i    轨道倾角，单位：弧度 

Ω   升交点赤经，单位：弧度 

ω   近地点幅角，单位：弧度 

f    真近点角(或M平近点角或E偏近点角)，单位：弧度 

注：角度转换为弧度的方式为：1角度的弧度＝1角度/180度乘以π 

见图14。有了这些参数，就能在地心赤道惯性坐标系(原点O在地心，X轴指向春分点方向，Z轴指向北极，OXYZ构成右手坐标系)下表示出二体卫星的位置和速度随真近点角的变化规律。公式如下： 

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\frac{a(1-e^2)}{1+e\cos f}\left[\begin{array}{c}\cos \mathrm{\Ω}\cos (\mathrm{\ω}+f)-\sin \mathrm{\Ω}\sin (\mathrm{\ω}+f)\cos i\\ \sin \mathrm{\Ω}\cos (\mathrm{\ω}+f)+\cos \mathrm{\Ω}\sin (\mathrm{\ω}-f)\cos i\\ \sin (\mathrm{\ω}+f)\sin i\end{array}\right]$

其中将如下三角函数展开： 

cos(ω+f)＝cosωcosf-sinωsinf 

sin(ω+f)＝sinωcosf-cosωsinf 

$\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]=$

$\sqrt{\frac{\mathrm{\μ}}{a(1-e^2)}}(-\sin f\left[\begin{array}{c}\cos \mathrm{\ω}\cos \mathrm{\Ω}-\sin \mathrm{\ω}\sin \mathrm{\Ω}\cos i\\ \cos \mathrm{\ω}\sin \mathrm{\Ω}+\sin \mathrm{\ω}\cos \mathrm{\Ω}\cos i\\ \sin \mathrm{\ω}\sin i\end{array}\right]+$

$(e+\cos f)\left[\begin{array}{c}-\sin \mathrm{\ω}\cos \mathrm{\Ω}-\cos \mathrm{\ω}\sin \mathrm{\Ω}\cos i\\ -\sin \mathrm{\ω}\sin \mathrm{\Ω}+\cos \mathrm{\ω}\cos \mathrm{\Ω}\cos i\\ \cos \mathrm{\ω}\sin i\end{array}\right])$

在此基础上，如果给出卫星经过近地点的时刻，通过求解开普勒方程，就能给出卫星随着时间的变化规律。为了得到真近点角和时间的关系，还需要引入两个随时间有关的角度量： 

M    平近点角 

E    偏近点角 

其中 

M＝n(t-t_p) 

它指的是卫星一个周期内的某个时刻经过的平均角度，t是当前时刻，tp是卫星过近地点的时刻，n是卫星的平均轨道角速度，它的公式是： 

而E是偏近点角，它和真近点角F的关系如图15所示： 

数学表达式为： 

$\tan \left(\frac{f}{2}\right)=\sqrt{\frac{1+e}{1-e}}\tan \left(\frac{E}{2}\right)$

而偏近点角E和平近点角M之间的关系可以用开普勒方程来表示： 

M＝n(t-t_p)＝E-esinE 

由上面的分析可以知道，通过求解具有非线性的开普勒方程，得到E，再通过E和f之间的关系，就可以得到真近点角f和时间t之间的对应关系。因为开普勒方程的非线性特性，因此只能进行非线性的求解，其中的一些方法如下： 

①.一阶微分法 

通过求微分，可以得到真近点角的变化步长和时间的变化步长之间的关系如下： 

$\sqrt{\frac{\mathrm{\μ}}{a^3}}\mathrm{\Δt}=\frac{(1-e\frac{e+\cos f}{1+e\cos f})\sqrt{1-e^2}}{1+e\cos f}\mathrm{\Δf}$

其中μ＝398600.44km³/s²(注意单位的一致性问题) 

②.级数计算式 

上述公式中的方法来进行时间的积分，误差较大，可改为采用如下方法： 

$f=M+(2e-\frac{e^3}{4})\sin M+\frac{5}{4}{e}^2\sin 2M+\frac{13}{12}{e}^3\sin 3M$

其中M＝n(t-t₀) 

③.牛顿迭代法 

当采用上述方法还不够精确的时候，可采用牛顿法进行迭代求解，求解方法如图16： 

f(E)＝(E-esinE)-M 

$\frac{\mathrm{df}\left(E_k\right)}{\mathrm{dt}}=1-e\cos E_k$

则有： 

$E_{k+1}=E_k-\frac{f\left(E_k\right)}{f\left(E_k\right)/\mathrm{dt}}$

基于以上公式就可以生成一条二体轨道曲线，其流程如图17所示： 

以上模块可以简计为： 

由此可以编制卫星类satellite_calculate中的TwoBody()函数，运行该函数并在平台上显示可得： 

使用构造函数satellite_calculate(char s[20])调用input.txt文件，文件中内容为： 

7000 0.01 30 20 60 0 0  24 60 279.306 

格式为： 

satellite_calculate satellite1(″satellite1_TwoBody″)； 

然后调用成员函数TwoBody()可得 

输出文件夹为satellite1_TwoBody的文本数据集合，在平台上作图，可以得到如下的低轨卫星轨道如图18所示： 

其星下点轨迹为图19所示： 

使用构造函数satellite_calculate(char s[20]，double da，double de，double di，double d_RAAN，double dw，double dt0，double dtp，double dtf，double dData_t，double dG0)，设置参数为： 

satellite_calculate satellite2(″satellite2_TwoBody″，500，20000，70，20，60，0，0，72，60，279.306)； 

则输出文件夹为satellite1_TwoBody的文本数据集合，在平台上作图，可得如下高轨卫星轨道：如图20所示。其星下点轨迹为图21所示： 

卫星的轨道周期可以表示为： 

$T=2\mathrm{\π}\sqrt{\frac{a^3}{\mathrm{\μ}}}$

由此可以看出，卫星的轨道周期和轨道半长轴有关系，调整半长轴就能设计不同周期的卫星轨道。为了使卫星在重复了若干圈之后其行下点轨迹再一次重合于原来的轨迹，只需要适当调整半长轴使轨道周期的整数倍和地球自转周期的整数倍相等即可。由此而产生的轨道称之为周期轨道。如图22所示： 

$T=\frac{24}{13}\×60\×60=6646.153s$

则有： 

$a=\sqrt[3]{\frac{T^2\mathrm{\μ}}{4{\mathrm{\π}}^2}}=7640.23\mathrm{km}$

计算可得： 

④卫星的地心赤道惯性坐标系下位置坐标在地理坐标系下的表示 

选择地球标准椭球体，知道地球的扁率和赤道半径a_e，便可以求出地球的极半径b_e，三者的相互关系如下： 

$f=\frac{a_e-b_e}{a_e}$

在此基础上在大地直角坐标系当中，可以得到地球椭球体的模型为： 

$\frac{{x1}^2}{a_e^2}+\frac{{y1}^2}{a_e^2}+\frac{{z1}^2}{b_e^2}=1$

所有地球椭球体表面上的点都将满足这个方程，又因为在知道了经纬(L)，纬度(B)，高程H的情况下： 

$B=\arctan \left(\frac{\tan \mathrm{\δ}}{{(1-f)}^2}\right)$

x1＝r_eBcosδcosL 

y1＝r_eBcosδsinL 

z1＝r_eBsinδ 

r＝r_eB+H 

整理可得： 

$r_{\mathrm{eB}}=\sqrt{\frac{1}{\frac{1}{a_e^2}+(\frac{1}{b_e^2}-\frac{1}{a_e^2}){\sin }^2\mathrm{\δ}}}$

又设卫星在惯性坐标系下的位置为(x，y，z)，因为： 

$\sin \mathrm{\δ}=\frac{z}{\sqrt{x^2+y^2+z^2}}$

所以卫星在某位置时其星下点处的地球半径： 

$r_{\mathrm{e\δ}}=\sqrt{\frac{a_e^2{b}_e^2(x^2+y^2+z^2)}{b_e^2(x^2+y^2+z^2)+(a_e^2-b_e^2)z^2}}$

则有高程为： 

$H=\sqrt{(x^2+y^2+z^2)}-r_{\mathrm{e\δ}}$

又卫星位置的赤经和赤纬可以表示为： 

$\mathrm{\α}=\arctan \left(\frac{y}{x}\right)$

$\mathrm{\δ}=\arcsin \left(\frac{z}{\sqrt{x^2+y^2+z^2}}\right)$

则经度可以表示为： 

L＝α-[G₀+ω_e(t-t₀)] 

地理纬度可以表示为： 

$B=\arctan \left(\frac{\tan \mathrm{\δ}}{{(1-f)}^2}\right)$

⑤从地理坐标系当中调取数据到地心赤道惯性坐标系 

α＝L+[G₀+ω_e(t-t₀)] 

tanδ＝(1-f)²tanB 

$\cos \mathrm{\δ}=\sqrt{\frac{1}{1+{\tan }^2\mathrm{\δ}}}$

$\sin \mathrm{\&delta;}=\left\{\begin{array}{cc}\sqrt{1-{\cos }^2\mathrm{\&delta;}}& \mathrm{\&delta;}>0\\ -\sqrt{1-{\cos }^2\mathrm{\&delta;}}& \mathrm{\&delta;}<0\end{array}\right.$

$r_{\mathrm{e\&delta;}}=\sqrt{\frac{1}{\frac{1}{a^2}+(\frac{1}{b^2}-\frac{1}{a^2}){\sin }^2\mathrm{\&delta;}}}$

r＝r_eδ+H 

则有： 

$\left\{\begin{array}{c}x=r\cos \mathrm{\&delta;}\cos \mathrm{\&alpha;}\\ y=r\cos \mathrm{\&delta;}\sin \mathrm{\&alpha;}\\ z=r\sin \mathrm{\&delta;}\end{array}\right.$

⑥卫星位置速度方式到轨道根数方式的转化 

卫星的位置和速度可以由以下公式来表示： 

$\stackrel{\&RightArrow;}{r}=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$ $\stackrel{\&RightArrow;}{v}=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]$

则设轨道动量矩为 $h=\left[\begin{array}{c}{h}_x\\ h_y\\ h_z\end{array}\right],$ 则有： 

$\left\{\begin{array}{c}{h}_x={\mathrm{yv}}_z-{\mathrm{zv}}_y\\ h_y={\mathrm{zv}}_x-{\mathrm{xv}}_z\\ h_z={\mathrm{xv}}_y-{\mathrm{yv}}_x\end{array}\right.$

设轨道单位质量能量为E，则有 

$E=\frac{V^2}{2}-\frac{\mathrm{\&mu;}}{r}$

其中 $r=\sqrt{x^2+y^2+z^2}$ $v=\sqrt{v_x^2+v_y^2+v_z^2}$

节线N为： 

$\stackrel{\&RightArrow;}{N}=\left\{\begin{array}{c}-\frac{h_y}{\sqrt{h_x^2+h_y^2}}\\ \frac{h_x}{\sqrt{h_x^2+h_y^2}}\\ 0\end{array}\right\}$

偏心率矢量e可以表示为： 

$\stackrel{\&RightArrow;}{e}=\left[\begin{array}{c}{e}_x\\ e_y\\ e_z\end{array}\right]=\left[\begin{array}{c}\frac{v_y{h}_z-v_z{h}_y}{\mathrm{\&mu;}}-\frac{x}{r}\\ \frac{v_z{h}_x-v_x{h}_z}{\mathrm{\&mu;}}-\frac{y}{r}\\ \frac{v_x{h}_y-v_y{h}_x}{\mathrm{\&mu;}}-\frac{z}{r}\end{array}\right]$

则有轨道倾角为： 

$i=\arccos \frac{h_z}{h}$

升交点赤经为： 

$\tan \mathrm{\&Omega;}=\frac{h_x}{-h_y}$

$\cos \mathrm{\&Omega;}=\frac{-h_y}{\sqrt{h_x^2+h_y^2}}$

当sinΩ＞0时： 

$\mathrm{\&Omega;}=\mathrm{ar}\cos \frac{-h_y}{\sqrt{h_x^2+h_y^2}}$

当sinΩ＜0时： 

$\mathrm{\&Omega;}=2\mathrm{\&pi;}-\mathrm{ar}\cos \frac{-h_y}{\sqrt{h_x^2+h_y^2}}$

轨道半长轴为： 

$a=-\frac{\mathrm{\&mu;}}{2E}$

偏心率为： 

$e=\sqrt{1-\frac{h^2}{\mathrm{\&mu;a}}}$

真近点角为： 

(ye_z-ze_y)h_x+(ze_x-xe_z)h_y+(xe_y-ye_x)h_z＜0时： 

$f=\arccos \left(\frac{h^2/\mathrm{\&mu;}-r}{\mathrm{re}}\right)$

(ye_z-ze_y)h_x+(ze_x-xe_z)h_y+(xe_y-ye_x)h_z＞0时： 

$f=2\mathrm{\&pi;}-\arccos \left(\frac{h^2/\mathrm{\&mu;}-r}{\mathrm{re}}\right)$

近地点幅角： 

当 $-{\mathrm{zh}}_x^2-{\mathrm{zh}}_y^2+({\mathrm{xh}}_x+{\mathrm{yh}}_y)h_z<0$ 时： 

$\mathrm{\&omega;}=\arccos (\frac{-{\mathrm{xh}}_y}{r\sqrt{h_x^2+h_y^2}}+\frac{{\mathrm{yh}}_x}{r\sqrt{h_x^2+h_y^2}})-f$

当 $-{\mathrm{zh}}_x^2-{\mathrm{zh}}_y^2+({\mathrm{xh}}_x+{\mathrm{yh}}_y)h_z>0$ 时： 

$\mathrm{\&omega;}=2\mathrm{\&pi;}-\arccos (\frac{-{\mathrm{xh}}_y}{r\sqrt{h_x^2+h_y^2}}+\frac{{\mathrm{yh}}_x}{r\sqrt{h_x^2+h_y^2}})-f$

⑥无摄二体轨道动力学数值积分方法求解初始位置速度问题 

在无摄无控条件下，轨道动力学可以由如下方程组表示： 

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\frac{\mathrm{\&mu;x}}{r^3}\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\frac{\mathrm{\&mu;y}}{r^3}\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\frac{\mathrm{\&mu;z}}{r^3}\end{array}\right.$

假如设： 

$y=\left[\begin{array}{c}x\\ y\\ z\\ v_x\\ v_y\\ v_z\end{array}\right]$

$f(t,y)=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -\frac{\mathrm{\&mu;y}}{r^3}\\ -\frac{\mathrm{\&mu;y}}{r^3}\\ -\frac{\mathrm{\&mu;z}}{r^3}\end{array}\right]$

常微分方程的求解方法： 

欧拉方法： 

y_n+1＝y_n+hf(t_n，y_n) 

四阶龙格库塔方法： 

$y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)$

K₁＝f(t_n，y_n) 

$K_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}{K}_1)$

$K_3=f(t_n+\frac{h}{2},y_n+\frac{h}{2}{K}_2)$

K₄＝f(t_n+h，y_n+hK₃) 

⑥摄动问题 

1)、摄动方程组的表示 

$\left\{\begin{array}{c}\frac{\mathrm{da}}{\mathrm{dt}}=\frac{2}{n\sqrt{1-e^2}}[F_{r}e\sin f+F_t(1+e\cos f)]\\ \frac{\mathrm{de}}{\mathrm{dt}}=\frac{\sqrt{1-e^2}}{\mathrm{na}}[F_r\sin f+F_t(\cos E+\cos f)]\\ \frac{\mathrm{d\&Omega;}}{\mathrm{dt}}=\frac{r\sin (\mathrm{\&omega;}+f)}{{\mathrm{na}}^2\sqrt{1-e^2}\sin i}{F}_n\\ \frac{\mathrm{di}}{\mathrm{dt}}=\frac{r\cos (\mathrm{\&omega;}+f)}{{\mathrm{na}}^2\sqrt{1-e^2}}{F}_n\\ \frac{\mathrm{d\&omega;}}{\mathrm{dt}}=\frac{\sqrt{1-e^2}}{\mathrm{nae}}[-F_r\cos f+F_t\frac{2+e\cos f}{1+e\cos f}\sin f]-\cos i\frac{r\sin (\mathrm{\&omega;}+f)}{{\mathrm{na}}^2\sqrt{1-e^2}\sin i}{F}_n\\ \frac{\mathrm{dM}}{\mathrm{dt}}=n-\frac{1-e^2}{\mathrm{nae}}[F_r(\frac{2\mathrm{er}}{p}-\cos f)+F_t(1+\frac{r}{p})\sin f]\end{array}\right.$

其中p为半通径，公式为： 

p＝a(1-e²) 

其中r为半径，公式为： 

$r=\frac{a(1-e^2)}{1+e\cos f}$

其中F_r，F_n，F_n为摄动力，考虑不同的摄动因素，摄动力的形式不相同 

2)、近地轨道地球形状摄动 

J₂摄动： 

$\left\{\begin{array}{c}{F}_r=-\frac{3}{2}{J}_2\frac{\mathrm{\&mu;}{R}_e^2}{r^4}[1-3{\sin }^2{i\sin }^2(\mathrm{\&omega;}+f)]\\ F_t=-\frac{3}{2}{J}_2\frac{\mathrm{\&mu;}{R}_e^2}{r^4}{\sin }^{2}i\sin 2(\mathrm{\&omega;}+f)\\ F_n=-\frac{3}{2}{J}_2\frac{\mathrm{\&mu;}{R}_e^2}{r^4}\sin 2i\sin (\mathrm{\&omega;}+f)\end{array}\right.$

对于J2摄动问题，如果仅仅考虑积分以后会按周期性进行积累误差的长周期摄动，，则摄动方程组可以简化为如下形式： 

$\left\{\begin{array}{c}\frac{\mathrm{da}}{\mathrm{dt}}=0\\ \frac{\mathrm{de}}{\mathrm{dt}}=0\\ \frac{\mathrm{di}}{\mathrm{dt}}=0\\ \frac{\mathrm{d\&Omega;}}{\mathrm{dt}}=-\frac{3}{2}\frac{J_2{a}_E^2}{p^2}n\cos i\\ \frac{\mathrm{d\&omega;}}{\mathrm{dt}}=\frac{3}{2}\frac{J_2{a}_E^2}{p^2}n(2-\frac{5}{2}{\sin }^{2}i)\\ \frac{\mathrm{dM}}{\mathrm{dt}}=n-\frac{3}{2}\frac{J_2{a}_E^2}{p^2}n(1-\frac{3}{2}{\sin }^{2}i)\sqrt{1-e^2}\end{array}\right.$

在Satellite_calculate类中编写成员函数J2_Perturbation()，使用构造函数satellite_calculate(char s[20])调用input.txt文件，文件中内容为： 

7000 0.01 30 20 60 0 0 60 60 279.306 

格式为： 

satellite_calculate satellite1(″satellite1_J2Perturbation″)； 

然后调用J2_Perturbation()函数，可得输出文件夹为satellite1_J2Perturbation的文本数据集合，在平台上调用数据可得：如图23 

一个月以后的轨道为图24。 

其他以此类推，最后通过实际观测数据来校准算法。 

Claims (9)

1.一种基于国产化技术开发的

卫星仿真开发与运行平台系统，其中，

该系统分别具有：

卫星轨道分析与三维空间任务显示、覆盖分析、导弹弹道预测分析、拦截器飞行工具分析、接近碰撞分析、来袭目标(弹道导弹)分析、链路分析、通信分析、电磁态势分析、雷达分析、推演分析等组成的卫星仿真开发与运行平台功能。

2.根据权利要求1所述的国产化技术开发，其中，

卫星仿真开发与运行平台(Virtual PrimePlatformSatellite Tool Kit)是中国研制生产的、底层用C/C++开发的、纯国产化的、基于国产组件化开发框架结构的+卫星仿真中间件组件、集成的一套卫星仿真开发包，打破了美国对华封锁禁售的美国Analytical Graphics公司的STK(Satellite ToolKit)卫星仿真开发包，在中国Satellite Tool Kit(STK)还是没有替代品。

3.根据权利要求1所述国产化技术，其中，

解决了国外同类产品(美国Analytical Graphics公司的Satellite Tool Kit)卫星仿真软件不支持中文字符显示的难题，以及支持国产二维GIS和三维GIS的难题。

4.根据权利要求1所述的国产化技术，其中，

底层采用C/C++语言开发，解决了跨操作系统平台的难题，可在国产linux操作系统和国产嵌入式实时操作系统上运行。

5.根据权利要求2所述的国产

组件化开发框架结构，其中

采用先进的开放式框架，实现灵活的软件框架结构：内核采用主控/组件与插件相结合的开放式总体架构；每个组件完成一个完整的业务，通过添加组件进行业务的扩展；每个插件点完成一个业务环节(一类算法)，通过增加插件进行算法的扩展；采用主控+注册组件+动态库形式，实现软件结构可伸缩的技术架构，软件功能可随时扩展、大小结构可伸缩，便于不同权限级别用户层次划分部署；

基于国产

组件化开发框架结构、采用组件化开发模式，便于不同单位、部门、项目组之间的协同开发与集成加载，采用多层次主控/组件结构，由主控(只有一个)、二级主控(可能多个)及若干组件组成，控主要完成组件注册、消息分发、组件通信、接口查询等的功能，组件完成具体的业务，可以接收并处理主控传来的各种命令、消息，并提供功能接口并外部调用，二级主控既是主控的一个组件，又是其管理的组件的主控。

6.根据权利要求3所述的支持国产二维GIS和三维GIS的难题，其中，

与国军标无缝集成：可集成国产Eagle VVP-DTIS二维GIS显示引擎组件v1.0、可使用总参测绘局GIS数据包括河流、省区界、居民区、铁路、机场、道路、高程等信息；

与国产系统无缝集成：可集成国产Eagle VVP-3D三维视景仿真引擎组件集v1.0、可加载Deep Eye三维战场态势显示平台v3.0；

与第三方软件接口：能够与MATLAB、STK、HLA、DDS、VegaPrime、WorldWind做接口或集成，同时能够与总参二维GIS和军用一体化平台三维GIS无缝集成。

7.根据权利要求4所述的底层采用C/C++语言开发，其中，

支持跨操作系统平台：提供windows XP/windows 7、VxWorks5.5/6.6、Linux、国产红旗Linux操作系统的跨平台编译支持，

底层平台与应用无关，中间件层由专业模型组件和功能组件组成，应用层与用户业务定制相关，实现模型性能验证框架、模型组件、模型参数、模型数据(文件)四者相分离。

8.根据权利要求1所述的卫星轨道分析与三维空间任务显示，其中，

支持空间任务的多视角观测：支持卫星观测包括对地观测(扫描地球)、恒星观测、月球观测、冷空间观测等，能够配合数据驱动模块，产生生动逼真的三维空间探测任务显示。

9.根据权利要求1所述的卫星轨道分析与三维空间任务显示、覆盖分析、导弹弹道预测分析、拦截器飞行工具分析、接近碰撞分析、来袭目标(弹道导弹)分析、链路分析、通信分析、电磁态势分析、雷达分析、推演分析，其中，这些都是美国对华封锁禁售的模块分析工具。

Images