卫星组网区域覆盖效能的仿真及计算方法与流程

本发明涉及卫星仿真，具体涉及一种卫星组网区域覆盖效能的仿真及计算方法。

背景技术：

1、大规模卫星组网对地面区域的时空覆盖效能仿真是卫星组网化、体系化应用效能仿真的基础，由于这类仿真涉及的卫星数量众多、任务区域范围大，导致传统的串行仿真评估算法计算量非常大，无法满足高性能仿真运行、动态快速指标评估等需求。

技术实现思路

1、有鉴于此，本发明旨在提出一种卫星组网区域覆盖效能的仿真及计算方法，解决现有仿真系统中大规模、多任务卫星组网对大面积地面区域的遥感、通信、导航等覆盖效能计算效率不高的问题，在保证计算精度的同时有效减少计算时间。

2、本发明实施例提供一种卫星组网区域覆盖效能的仿真及计算方法，所述方法包括：

3、s100，执行仿真场景编辑与初始化；

4、s200，通过多个仿真cpu线程，并行仿真卫星组网的每个卫星的空间状态，其中，卫星组网的数量为一个或多个；

5、s300，通过至少一个计算cpu线程，将步骤s100中执行仿真场景编辑与初始化后获得的数据，以及步骤s200中获得的每个卫星的空间状态数据，写入gpu；

6、s400，通过gpu并行计算每个卫星组网的覆盖效能指标。

7、优选地，步骤s100包括：

8、设置仿真持续时间、仿真步长、场景初始纪元时间和卫星初始状态；

9、选择卫星组网任务类型，卫星组网任务类型包括遥感、通信和导航三种；

10、划分目标区域，目标区域的数量至少为一个，每个卫星组网与至少一个目标区域具有对应关系；以及

11、将每个卫星组网与其对应的目标区域关联；

12、其中，所述划分目标区域包括：设置目标区域的网格划分的精度值，每个目标区域按照设置的精度值进行网格点划分，每个网格点由经纬度表示，并通过每个目标区域的多边形顶点的经纬度集合来描述相应目标区域的边界范围。

13、优选地，在步骤s200中，每个卫星均对应一个仿真cpu线程。

14、优选地，在步骤s200中，仿真单个卫星的空间状态，包括：

15、计算卫星的瞬时轨道六根数、惯性系位置、星下点位置、卫星轨道高度和卫星本体姿态；

16、当卫星的任务类型为遥感类任务时，还需要计算上一步仿真时刻tk-1到当前仿真时刻tk的卫星传感器视场对地覆盖范围，卫星传感器视场对地覆盖范围为{l(tk-1),l(tk),r(tk),r(tk-1)}四个顶点构成的矩形，每个顶点由地球固连坐标系下的经纬度表示，其中，l(tk)和r(tk)为卫星传感器在tk时刻横向视场边界视线与地球的两个交点，l(tk-1)和r(tk-1)为卫星传感器在tk-1时刻横向视场边界视线与地球的两个交点。

17、优选地，步骤s300包括：

18、s310，根据卫星组网的数量，调用相同数量的计算cpu线程，每个计算cpu线程对应一个卫星组网，且每个卫星组网对应一个目标区域；

19、s320，多个计算cpu线程并行运行，排队调用gpu设备，每个计算cpu线程均将对应的卫星组网的数据写入相应gpu；

20、写入gpu的数据包括：

21、目标区域的所有网格点的经纬度信息；

22、卫星组网每个卫星的惯性系位置、星下点位置和卫星轨道高度；

23、对于遥感类任务的卫星，还包括卫星传感器视场对地覆盖范围；

24、对于导航类和通信类任务的卫星，还包括最小观测仰角。

25、优选地，在步骤s400中，计算每个卫星组网的覆盖效能指标包括：

26、s410，基于cuda进行第一次多线程网格划分，所有gpu线程并行调用第一个内核函数，对tk时刻单个卫星sati(i∈[1,m])对于目标区域单个网格点pj(j∈[1,n])的可见性指标v(sati,pj,tk)进行计算，将计算得到的每个卫星对于每个网格点的可见性指标集合写入gpu的全局内存；

27、s420，根据步骤s410中第一个内核函数的计算结果，基于cuda进行第二次多线程网格划分，所有gpu线程并行调用第二个内核函数，对卫星组网对于目标区域每个网格点的覆盖效能指标进行计算；

28、其中，对于导航类任务的卫星组网，还需要计算卫星组网对于目标区域每个网格点的导航精度因子；

29、s430，根据步骤s410和步骤s420中两个内核函数的计算结果，调用gpu的第三个内核函数，采用归约算法，对卫星组网对目标区域{p1,p2,...,pj,...,pn}(j∈[1,n])的整体覆盖效能指标进行统计计算；

30、其中，对于导航类任务的卫星组网，还需要计算卫星组网对于目标区域的整体导航精度因子指标；

31、s440，gpu并行处理结束，将所有计算结果写回cpu。

32、优选地，在步骤s410中，包括：

33、对于遥感类任务的卫星，可见性指标v(sati,pj,tk)的计算方式为：判别点pj是否在卫星传感器sati视场覆盖{l(tk-1),l(tk),r(tk),r(tk-1)}多边形区域内，判别方法采用水平/垂直交叉点数判别法；

34、对于通信类和导航类任务的卫星，可见性指标v(sati,pj,tk)的计算方式为：设地面网格点的经纬度分别为λ0、φ0，卫星星下点的经纬度分别为λs、φs，观测点对卫星的最小观测仰角为e，则星下覆盖区对应的地球中心角γ为：

35、

36、仰角e为：

37、

38、其中，re为地球半径，h为卫星轨道高度，emin为卫星最小观测仰角，则可见性指标v(sati,pj,tk)函数为：

39、

40、优选地，在步骤s420中，包括：

41、计算在tk时刻，卫星组网{sat1,sat2,...,sati,...,satm}(i∈[1,m])整体对目标区域{p1,p2,...,pj,...,pn}中每个网格点的瞬时可见性和瞬时覆盖重数两个指标；

42、其中，在tk时刻，卫星组网对单个网格点pj(j∈[1,n])的瞬时可见性指标表示为c(pj,tk)，瞬时覆盖重数指标表示为nfc(pj,tk)，计算结束后将每个网格点的可见性指标和瞬时覆盖重数指标写入gpu全局内存，计算公式分别为：

43、

44、

45、统计在(0,tk]时间段内，卫星组网{sat1,sat2,...,sati,...,satm}(i∈[1,m])整体对目标区域{p1,p2,...,pj,...,pn}中每个网格点pj(j∈[1,n])的可见时刻集tws(pj,tk)，统计结束后将每个网格点的可见时刻集写入gpu全局内存，tws(pj,tk)可表示为：

46、

47、其中，tw(pj,tp,tq)表示在(0,tk]时间段内可见性连续为1的一个时间段子集，可表示为：

48、tw(pj,tp,tq)＝{(tp,tq]|c(pj,tp+1)×c(pj,tp+2)×……×c(pj,tq-1)×c(pj,tq)＝1}；

49、则在(0,tk]时间段内卫星组网对点pj的覆盖次数nc(pj,tk)等于可见时刻集tws(pj,tk)的元素个数，其公式为：

50、nc(pj,tk)＝|tws(pj,tk)|；

51、基于仿真推演，已知tk-1时刻卫星组网对点pj的可见时刻集表示为：

52、

53、则tk时刻卫星组网对点pj的可见时刻集按如下公式进行计算：

54、

55、统计在(0,tk]时间段内，卫星组网{sat1,sat2,...,sati,...,satm}(i∈[1,m])整体对目标区域{p1,p2,...,pj,...,pn}(j∈[1,n])每个网格点的累计覆盖效能指标，计算结束后将所有网格点的相应效能指标写入gpu全局内存，具体指标包括：总覆盖时间、覆盖时间百分比、累计平均覆盖时间、累计最大覆盖时间、累计最大覆盖间隙、累计平均覆盖间隙、累计最大重访时间和累计最小/最大覆盖重数；

56、其中：

57、总覆盖时间是指在(0,tk]仿真时间段内，卫星组网对地面网格点可见性为1的持续时间长度总和，设在(0,tk]时间段内卫星组网对点pj的总覆盖时间表示为ttc(pj,tk)，则ttc(pj,tk)的计算公式是：

58、

59、由于每一步的仿真计算基于上一步的仿真计算结果，因此ttc(pj,tk)每步仿真的计算公式为：

60、ttc(pj,tk)＝ttc(pj,tk-1)+(tk-tk-1)×c(pj,tk)；

61、覆盖时间百分比是指在(0,tk]仿真时间段内，卫星组网对地面网格点的覆盖时间百分比等于其总覆盖时间除以时间总长tk，设(0,tk]时间段内卫星组网对点pj的覆盖时间百分比表示为tpc(pj,tk)，则tpc(pj,tk)的计算公式是：

62、

63、累计平均覆盖时间的计算方式为：设(0,tk]时间段内卫星组网对点pj的平均覆盖时间表示为tmeanc(pj,tk)，则累计平均覆盖时间为：

64、

65、累计最大覆盖时间等于(0,tk]时间段内最长的连续覆盖时间，卫星组网对点pj的最大覆盖时间表示为tmaxc(pj,tk)，则累计最大覆盖时间计算公式是：

66、

67、累计最大覆盖间隙表示卫星组网对地面点连续不可见的最长时间，设在(0,tk]时间段内卫星组网对点pj的最大覆盖间隙表示为tmaxcg(pj,tk)，则累计最大覆盖间隙计算公式是：

68、

69、其中m＝nc(pj,tk)，是可见时刻集tws(pj,tk)按时间从小到大排序的第i个连续可见时间段子集；

70、累计平均覆盖间隙表示卫星组网对地面点覆盖中断的平均时间长度，等于覆盖中断总时间除以覆盖中断次数，设在(0,tk]时间段内卫星组网对点pj的最大覆盖间隙表示为tmeancg(pj,tk)，则累计平均覆盖间隙计算公式是：

71、

72、其中m＝nc(pj,tk)，是可见时刻集tws(pj,tk)按时间从小到大排序的第i个连续可见时间段子集；

73、累计最大重访时间是卫星对地面点的重访时间是指从上一次连续覆盖的起始时刻到这一次连续覆盖的起始时刻的时间间隔，最大重访时间是指这个时间间隔的最大值，在(0,tk]时间段内卫星组网对点pj的最大重访时间tmaxrt(pj,tk)的计算公式为：

74、

75、其中m＝nc(pj,tk)，是可见时刻集tws(pj,tk)按时间从小到大排序的第i个连续可见时间段子集，i＝1,2,……,nc(pj,tk)且1≤pi＜qi≤k；

76、累计最小/最大覆盖重数是统计在(0,tk]时间段内卫星组网对点pj的每个时刻覆盖重数的最小值、最大值和平均值，网格点pj的最小、最大覆盖重数分别表示为min_nfc(pj,tk)、max_nfc(pj,tk)，计算公式分别是：

77、min_nfc(pj,tk)＝min{nfc(pj,tt)|t＝1,2,……,k}；

78、max_nfc(pj,tk)＝max{nfc(pj,tt)|t＝1,2,……,k}；

79、对于导航类任务的卫星组网，则还需要计算tk时刻卫星组网对目标区域每个网格点的导航精度因子指标，导航精度因子包括空间精度因子、钟差精度因子和几何精度因子，计算结束后将每个网格点的导航精度因子指标写入gpu全局内存；假设有n个导航卫星进行组网，具体如下：

80、根据网格点位置、卫星位置计算协因数阵q：

81、

82、

83、

84、

85、其中，(x0,y0,z0)为网格点位置坐标，(xi,yi,zi)为历元t时导航卫星i在协议地球坐标系中的坐标分量；

86、分别计算空间精度因子pdop、钟差精度因子tdop和几何精度因子gdop为：

87、

88、

89、

90、优选地，在步骤s430中，包括：

91、采用归约算法，统计卫星组网对目标区域网格点的瞬时最大、最小和平均覆盖重数nfc_max(tk)、nfc_min(tk)、nfc_mean(tk)，公式如下：

92、

93、

94、

95、根据nfc_max(tk)、nfc_min(tk)、nfc_mean(tk)的计算公式计算总覆盖时间、覆盖时间百分比、累计平均覆盖时间、累计最大覆盖时间、累计最大覆盖间隙、累计平均覆盖间隙、累计最大重访时间、空间精度因子、钟差精度因子和几何精度因子的最大、最小和平均值，计算结束后将计算结果写入gpu全局内存；

96、采用归约算法，计算tk时刻卫星组网对整个区域的可见时刻集tws(tk)，整个区域的可见时刻集为每个网格点的可见时刻集的并集，表示为：

97、

98、对于区域网格点集合{p1,p2,...,pj,...,pn}(j∈[1,n])，计算在(0,tk]时间段内卫星组网对目标区域的累计覆盖百分比，公式为：

99、

100、计算结束后将目标区域的累计覆盖百分比指标写入gpu全局内存。

101、优选地，在所述第一次多线程网格划分中，线程网格x维线程块数量为：

102、

103、线程网格y维线程块数量为：

104、griddim.y＝threadidx.y＝卫星或传感器数量；

105、一个线程块x维线程数量为：

106、blockdim.x＝一个线程束大小的整数倍；

107、一个线程块y维线程数量为：

108、blockdim.y＝1；

109、线程网格x维线程总数量为：

110、threadcountx＝griddim.x*blockdim.x；

111、线程网格y维线程总数量：

112、threadcounty＝卫星或传感器数量；

113、在所述第二次多线程网格划分中，线程网格x维线程块数量为：

114、线程网格y维线程块数量为：

115、griddim.y＝threadidx.y＝1；

116、一个线程块x维线程数量为：

117、blockdim.x＝一个线程束大小的整数倍；

118、一个线程块y维线程数量为：

119、blockdim.y＝1；

120、线程网格x维线程总数量为：

121、threadcountx＝griddim.x*blockdim.x；

122、线程网格y维线程总数量：

123、threadcounty＝1。

124、本发明的卫星组网区域覆盖效能的仿真及计算方法能够支持多cpu+多gpu的异构多处理器混合体系结构，通过多cpu线程、多gpu设备、多gpu线程三级并行机制，有效利用新型处理器、通信、同步等机制优化仿真资源，实现卫星组网覆盖效能仿真的高效并行运行。通过将cuda架构与卫星组网覆盖效能不同层次指标的计算和统计方法有机结合，充分利用cuda并行存储、通信和计算机制，实现分层并行的高效卫星组网覆盖效能并行仿真和计算方法，保证计算精度的同时有效减少计算时间。