CN102201008A - 一种基于gpu的快速星表检索方法 - Google Patents

Abstract

本发明提供了一种基于GPU的星表检索方法，其步骤包括：1)采用经纬度划分法将星表划分为星区存储；2)计算星敏感器圆形视场区域所覆盖的星区，该过程先计算视轴的球面坐标；然后计算圆形视场的经纬度边界，在分区星表中查找圆形视场所覆盖的星区；3)将星区内恒星传入GPU全局内存以恒星集方式存储，GPU恒星计算线程网格划分为若干个线程块，每个线程块负责处理一个恒星集合，每一个线程负责对一个恒星进行检索；4)各线程检索出落入圆形视场内的恒星，在圆形视场内检索出落入CCD成像面的恒星。本发明方法将各类视场下的检索时间控制在毫秒级别，使得在大视场、宽星等域下满足星图生成等各类仿真应用对恒星检索的实时要求。

Description

一种基于GPU的快速星表检索方法

技术领域

本发明属于计算机仿真领域，具体涉及一种具有快速实时特点的星表检索方法，该检索方法能够在给定大小的视场下，快速检索出星表中落入视场内的恒星。

背景技术

星敏感器是卫星姿态确定与控制***中重要的星上设备，主要采用快速的星图识别算法以确定卫星的姿态。为了测试星敏感器的识别算法性能，需要使用星模拟器，模拟任意时刻，任意飞行器姿态下的星敏感器视场内的星图。计算机模拟星敏感器等航天设备拍摄到的背景星图，首先根据探测视场光轴的赤经、赤纬，以给定的星敏感视场大小以及视场绕光轴的旋转角度确定视场范围，从星表提取落入视场内的恒星，最后将这些星点以二维图像形式显示出来。落人星敏感器视场内的观测星与模拟天球分布的所有观测星相比只占极少数，这就存在一个在星表快速检索落入星敏感器视场的观测星的问题。设计高效率的星表检索算法是提高星图模拟速度，实现星空环境实时仿真的关键因素。

星表包含了一定数量恒星的星编号、星等、对应某一历元时刻的星位置(天赤经、天赤纬)等数据，构成了对空间天球坐标系下恒星的描述。星表检索以恒星星表为基础，通过给定视场大小和视轴指向确定搜索范围，从星表中检索出落入视场内的恒星。为了加快星表检索速度，文献[1](参考Bone Jeffery W.on-orbit star processing using multi-startrackers[J].Proc.SPIE，1994，22(21)：6-8)提出一种球矩阵结构来存储星表，通过优化星表结构以提高恒星存取效率，但该方法在矩阵维数较大时，存取时间并不能显著减少。文献[2](参考饶才杰，房建成.一种星图模拟中观测星提取的方法[J].光学精密工程，2004，12(2)：129-134.)提出了经纬度圈法划分星表，用赤经圈和赤纬圈将天球分成不重叠的区域，每个区域独立存储，该方法边界清晰，计算简单，然而在处理大规模星表时容易产生过多的分区，影响恒星查询效率；文献[3](参考张广军，魏新国，江洁.一种改进的三角形星图识别方法[J].航空学报，2006，27(6)：1150-1154.)通过构造天球的“内接正方体”将天球球面均匀分成6块，再细分块为更多子块，通过子块索引号读取恒星，这种方法查询方便，但对恒星的块边界判断较复杂，当恒星数量较大时，边界计算耗时较高。为了缩小搜索范围，文献[4](参考胡宜宁，巩岩.动态星图显示算法的设计与实现[J].宇航学报，2008，29(3)：849-853.)针对某一指向的视场，取视场光轴所在分区临近的八个子区作为搜索区，检索落入视场内的恒星，该方法有效缩小检索空间，然而只能针对特定大小的视场，这是因为天球上赤经弧长随纬度增加而逐渐缩短，当视轴指向不同天区，不同大小的视场覆盖的子区数目并不相同。文献[5](参考蔡志武，韩春好，陈金平.一种高精度星载导航星库的构建方法[J].测绘科学技术学报，2006，23(1)：29-32.)通过计算球表面上视场的经纬度边界范围，进而查找该视场所关联的子区，但该方法采用膨胀因子计算经纬度边界范围，计算精度不高。此外，已有的视场内恒星检索过程均在CPU上的串行实现，在恒星数量巨大时，检索效率并不高。随着星敏感器等航天设备观测能力的提升，可观测星等域越来越大。目前在大视场、宽探测星等域条件下，如何提高星表检索速度已经成为航天仿真在星空领域的亟需解决的问题。

发明内容

针对现有方法的缺点，本发明提出一种基于GPU的星表检索方法。该方法采用经纬度划分法对星表进行星区划分；提出了一种视场所覆盖的星区计算方法，解决视场边界范围计算不准问题；然后采用NIVIDIA提出的CUDA并行架构，设计基于GPU平台的视场恒星查找并行算法；具有快速检索出星表的效果。

为了实现本发明的目的，采用的技术方案概述如下：

一种基于GPU的星表检索方法，其步骤包括：

1)采用经纬度划分法将星表划分为星区存储；

2)计算星敏感器圆形视场区域所覆盖的星区，该过程先计算视轴的球面坐标；然后计算圆形视场的经纬度边界，在分区星表中查找圆形视场所覆盖的星区；

3)将星区内恒星传入GPU全局内存以恒星集方式存储，GPU恒星计算线程网格划分为若干个线程块，每个线程块负责处理一个恒星集合，每一个线程负责对一个恒星进行检索；

4)各线程检索出落入圆形视场内的恒星，在圆形视场内检索出落入CCD成像面的恒星。

所述步骤1)经纬度划分法是按照一定的经纬度间隔，用赤经圈和赤纬圈将天球分成不重叠的区域，每个区域内的恒星独立存储，并按照各个区域所在经度和纬度范围进行编号。

进一步，对同一纬度区间内的星区按照经度递增的顺序合并存储到一个星表文件，并且建立索引。

进一步，所述索引记录了各星区所在经度范围以及星区在星表文件中的恒星首地址、偏移长度。

所述步骤2)通过视场在天球的圆形投影面与赤经线构造球面三角形，通过球面几何关系计算出经纬度边界。

进一步，所述构造球面三角形方法为：以O为视轴与天球面的交点，过O点的经线与视场区域交与P₁、P₂两点，视场区域与天球两条相切经线切于Q₁、Q₂两点；所述球面几何关系为：

$\left\{\begin{array}{c}\cos (90-{\mathrm{\α}}_1)=\cos \mathrm{\θ}\cos (90-{\mathrm{\α}}_2)\\ \sin \mathrm{\θ}=\sin (90-{\mathrm{\α}}_1)\sin \mathrm{\Δ\δ}\end{array}\right.$

θ为视场半张角，α₁、α₂分别为O和Q₂的纬度，Δδ为Q₂和O的经度差。

所述步骤4)判断恒星落入圆形视场所在区域的方法是：计算出该恒星到观测点的连线与视轴的夹角，比较该夹角与视场张角的大小关系，如果该夹角大于视场张角，该恒星位于视场外，否则，恒星落入圆形视场所在区域。

所述步骤4)判断恒星坐标落入CCD像平面上的方法是：该恒星坐标在像平面的坐标范围内。

(1)本发明的积极效果

与过去的星表检索方法相比，本发明有如下优点：1)采用经纬度划分法对星表进行分割，并引入合并存储方式存储各星表分区，使得星表存储管理方便，且提高了查询效率；2)通过引入球面三角形法求出视场在星表天球所覆盖的精确经纬度边界，并以此为基础，得出视场所覆盖的星区，以这些星区作为检索区域而不是整个星表，缩小了搜索范围；3)设计了基于GPU平台的星表检索并行检索算法。利用当前GPU强大的数值计算能力和并行处理能力，将对海量恒星的查找过程进行多线程并行加速，获得了数十倍的加速比；4)将各类视场下的检索时间控制在毫秒级别，使得在大视场、宽星等域下满足星图生成等各类仿真应用对恒星检索的实时要求。

附图说明

图1为天球中视场边界表示图；

图2为视场与CCD成像面位置关系图；

图3为恒星检索并行模式示意图；

图4为GPU星表检索算法示意图；

图5为不同GPU平台星表检索的加速效果比较示意图。

具体实施方式

下面结合附图和具体实施例对本发明作详细的说明。

由于近年来GPU数值计算能力快速提升，可进行数值计算的并行处理核心已可达300个。设计基于GPU上的并行计算对数据规模大且可并行化的处理显示出强大优势，已经在几何造型、分子模拟、图像处理等领域取得了很好的加速效果。对星表检索的分析显示，视场内恒星的检索过程具有较高的并行度，非常适合使用GPU进行加速。这是因为，首先在大视场，宽探测星等域下，视场内恒星查找是一个计算密集、数据量大的处理过程，其计算量占整个星表检索计算量的85％～90％，高效且快速地实现该查找过程是提高整体检索速度的关键；其次，查找视场恒星的每一次处理都是对视场覆盖星区内的恒星进行一次独立的计算和判断，且该计算不会影响到其他恒星。这种典型的数据海量性和问题本身具备的良好并行结构，使得在GPU上进行并行设计具有天然优势。

因此，本发明给出具体的实施例，说明本发明基于GPU并行架构下的星表快速检索方法。首先，采用经纬度划分法对星表进行星区划分；其次，提出了一种探测视场所覆盖的星区计算方法，解决视场边界范围计算不准问题；然后采用NIVIDIA提出的CUDA并行架构，设计基于GPU平台的视场恒星查找并行算法；最后通过仿真实验验证了本算法的加速效果。

(1)星表划分

星表划分是星表检索的基础，它以视场大小为输入，确定分区粒度，采用经纬度划分法，生成可快速查询的分区星表。该分区星表可进行多次检索。

经纬度划分法类似地球分带，按照一定的经纬度间隔，用赤经圈和赤纬圈将天球分成不重叠的区域，每个区域内的恒星独立存储，并按照各个区域所在经度和纬度范围进行编号。然而，由于每个星区内的恒星独立存储，当星表密度较大时，分区间隔较小导致存储文件数过多，不易管理并且影响查询效率。

为了减少存储文件数，采用合并存储方式，对同一纬度区间内的星区按照经度递增的顺序合并存储到一个星表文件，并且建立索引。索引记录了各星区所在经度范围以及星区在星表文件中的恒星首地址、偏移长度。这样，在读取任一星区的恒星数据时，可先根据其所在的纬度区间查找到其对应的存储文件，再根据星区所在经度区间在索引中查找到文件中的首地址和偏移长度的记录，便可读取该星区所有的恒星。

(2)视场覆盖星区计算

视场覆盖星区查找，在分区星表的基础上，计算圆形视场区域所覆盖的星区。该过程先计算视轴的球面坐标；计算光轴球面坐标过程：

恒星在USNOA1.0星标中使用的是赤经、赤纬表示的球面坐标系，而相机光轴的方向表示采用的是地心惯性坐标系(笛卡尔坐标)，因此要进行将地心惯性坐标系转换为球面坐标系以计算相机视场的覆盖范围。相机的光轴笛卡尔坐标(x，y，z)与球面坐标(Ra，De)的关系为：

$\left\{\begin{array}{c}x=\mathrm{SGR}\·\cos \left(\mathrm{De}\right)\·\cos \left(\mathrm{Ra}\right)\\ y=\mathrm{SGR}\·\cos \left(\mathrm{De}\right)\·\sin \left(\mathrm{Ra}\right)\\ z=\mathrm{SGR}\·\sin \left(\mathrm{De}\right)\end{array}\right.$

其中Ra为赤经，范围[0，2π)，De为赤纬，范围[-π/2，π/2]，SGR为天球半径，理论上是无穷远，实际可以取一个较大的值。

然后计算星敏感器圆形视场的经纬度边界，在分区星表中查找视场所覆盖的分区。

针对圆形视场在天球中的经纬度边界的计算，通过视场在天球的圆形投影面与赤经线构造球面三角形，通过球面几何关系计算出任意指向下视场的精确边界。

如错误！未找到引用源。所示，考虑观测点与地球的垂直距离和与恒星的中心距离相比极其微小，为了方便，可以认为观测中心为地心。

Oe为观测中心，O为视轴与天球面的交点，过O点的经线与视场区域交与P₁、P₂两点，视场区域与天球两条相切经线切于Q₁、Q₂两点，则易知，计算出P₁、P₂的赤纬度，Q₁、Q₂的赤经度，即可确定视场纬度、经度范围。首先，先计算出轴点O的经纬度，然后通过O点纬度加上或减去视场的一半张角就可算出纬度范围，而计算Q₁、Q₂两点的经纬度坐标可通过求解球面三角形PQ₁O、PQ₂O得到，以球面三角形PQ₂O为例，构造球面直角三角形PQ₂O，与球心Oe组成球心三面角Oe-PQ₂O(图1右)，根据求解直角球面三角形的聂比尔定则，可列方程为：

$\left\{\begin{array}{c}\cos (90-{\mathrm{\α}}_1)=\cos \mathrm{\θ}\cos (90-{\mathrm{\α}}_2)\\ \sin \mathrm{\θ}=\sin (90-{\mathrm{\α}}_1)\sin \mathrm{\Δ\δ}\end{array}\right.---\left(1\right)$

式中，θ为视场半张角，α₁、α₂分别为O和Q₂的纬度，Δδ为Q₂和O的经度差。已知θ和α₁，解此方程，求出α₂和Δδ，即得到Q₂的经纬度坐标，Q₁也可同样求出。

计算出探测视场的经纬度范围后，根据星表中各星区的经纬度范围，找出其所覆盖的星区号。这样，对于星表中包含星区总数为N，视场覆盖星区为n，则搜索范围缩小为原来的n /N。

(3)视场内恒星检索

计算出探测视场所覆盖的星区后，只需以这些星区为检索范围剔除未落入视场的恒星，从而检索出落在视场内的恒星。在空间探测器视场中，观测星最终落入圆形视场内的方形CCD成像面中，二者位置关系如图2所示。

由于判断星落入圆形视场比成像面视场内要简单很多，同时避免星区中每颗星都进行是否落入成像面的复杂计算，先从星区中检索出在圆形视场内的恒星，然后在圆形视场内检索出落入CCD成像面的恒星。该过程中大部分均是较复杂的数值计算，在星表检索中耗时最多，约占整个计算量的90％，，也是优化的关键环节。可以看到，查找过程是对一定数量的恒星执行相同的计算过程，是典型的单指令、多数据计算，适合采用GPU进行并行加速。因此，将视场内恒星查找过程转换到GPU上进行并行计算，并设计CUDA高性能架构下的并行检索算法。

CUDA是NVIDIA推出的一种高性能GPU计算架构，充分发挥了GPU多核并行处理和强大数值计算能力。在CUDA架构下，一个程序由运行在CPU(主机)端的控制程序和GPU(设备)端上的内核程序(kernel)组成。CPU负责处理逻辑性强的事务处理和串行计算，GPU则专注高度线程化的并行处理任务；执行计算时，通常由主机端将准备好的数据传输到GPU内存中，内核程序计算完成后在将结果从GPU内存中取回。

在GPU运行时，大量并发线程同时并行执行内核程序。而这些线程按照两个并行层次组织。线程(Thread)先被组织成线程块(Block)，同一个线程块内的线程可以通过GPU芯片上的共享内存和同步机制进行通信；这些线程块又被组织成可以执行一个内核的线程网格(Grid)。不同线程块之间只能通过所有线程均可访问的全局内存来共享数据。设计视场内恒星查找并行计算架构时，为了使得每个线程对应一个恒星查找计算，要建立线程块到恒星集合的对应关系。基于此，在恒星传入GPU全局内存时，存储的方式不再按照星区存储，而是根据以线程块的大小为单位分成若干个恒星集，如图3所示。待处理恒星划分为若干个恒星集合，在部署GPU线程时，将恒星计算线程网格划分为若干个线程块，每个线程块负责处理一个恒星集合，线程块内的每一个线程负责对一个恒星进行检索程序，建立起线程到恒星的并行关系。在执行过程中，线程间执行互不影响，独立执行。

CUDA下恒星检索算法以分区星表和视场覆盖星区查找为基础，对视场内恒星查找过程转换到GPU上进行多线程并行加速。图4给出了CUDA下完整的并行星表检索算法。首先由CPU主程序负责星表检索的串行计算，包括计算视轴的球面坐标和计算视场覆盖的星区。通过计算得到待搜索星区，将星区内恒星传入GPU全局内存以恒星集方式存储，传输结束后启动GPU上的内核检索程序(kernel)，由设置好的线程网格并行执行，线程网格包含的总线程数与恒星数目一致，网格包含的线程块数与恒星集合一致，每个线程块内包含线程数目为256个。在线程启动后，各线程首先从全局内存中读取对应恒星信息，执行圆形视场检索，判断该恒星是否落入圆形视场所在区域，判断方法是计算出该恒星到观测点的连线与视轴的夹角，比较该夹角与视场张角的大小关系，如果该夹角大于视场张角，该恒星位于视场外，则停止操作；否则计算该恒星投影到CCD像平面的坐标，执行像平面检索。判断该恒星坐标是否落入CCD像平面上，判断方法是该恒星坐标是否在像平面的坐标范围内。待计算完成后，该线程将计算出的像平面内的恒星位置信息传入到GPU全局内存中。所有并行线程运行结束后，像平面上恒星信息由GPU全局内存传送到CPU内存中，完成星表检索。

(4)检索效率比较

在仿真应用中常用的星表如FK系列星表、SAO星表、第谷星表和USNO星表等。几种星表包含的恒星数和分布密度并不相同，以USNO星表收录的恒星数最多，分布密度也最大。为了测试算法对大规模星表检索的性能，实验选用USNO-AV1.0作为基础星表，该星表共包含了全天球488,006,860颗星。实验将可检索星等范围选取为3～15等星，并从2.5度到20度范围内分别取四个递增值作为实验视场大小(如表1所示，显示了视场大小和不同计算平台下检索时间)，涵盖了星空仿真中各类应用的需求。

实验环境为CPU采用Intel Xeon X5450 3.0GHz，内存4GB；为了验证加速效果，使用两块图形处理器测试GPU检索程序，分别为Quadro FX4600和GeForce GTX480，前者包含96个处理核，可同时并发9216个线程，后者包含300个处理核，可同时并发30720个线程。

本发明共设计三组仿真实验以验证星表检索加速算法：第一组实验运行基于CPU平台的传统检索算法，第二组实验运行基于GPU(Quadro FX4600)平台的并行检索算法，第三组实验运行在GPU(GeForce GTX480)平台并行检索算法。针对四个观测视场，多次运行三组检索程序，计算出各组实验的星表检索平均耗时，如表1所示。从数据对比可以看到，当视场较小时，相比CPU平台，GPU加速优势还不是十分显著。这是由于在视场较小时，待处理恒星数目相对较少，并发线程数也较少，并没有充分利用GPU并发处理能力。当视场增大，恒星数量快速增加，GPU平台并行处理优势显著。这是因为在CPU平台需要对每颗恒星顺序处理，检索时间迅速增加。而GPU平台下，恒星检索过程由大量的并行线程执行，取得了很高的性能提升。当视场大小为20.0*20.0时，CPU需要33.565ms，而GPU并行检索只需要1.153ms，已可满足各类星图仿真的实时要求。表2给出了GPU并行加速星表检索算法所获得的加速比。随着视场增大，加速比快速增加，可以看到，在视场大小为20.0*20.0时，基于GPU(GTX480)平台的检索速度比CPU平台要快近60倍。图5给出了FX4600和GTX480平台的加速效果对比。图中曲线走势明显的反应了二者加速效果差距。这是由于GTX480平台相比FX4600可同时并发更多的线程数，并且在计算能力和内存访问速度也具有较大优势，因此在GTX480平台上的星表检索加速性能也更优异。

表1

表2

Claims (8)

1.一种基于GPU的星表检索方法，其步骤包括：

1)采用经纬度划分法将星表划分为星区存储；

2)计算星敏感器圆形视场区域所覆盖的星区，该过程先计算视轴的球面坐标；然后计算圆形视场的经纬度边界，在分区星表中查找圆形视场所覆盖的星区；

3)将星区内恒星传入GPU全局内存以恒星集方式存储，GPU恒星计算线程网格划分为若干个线程块，每个线程块负责处理一个恒星集合，每一个线程负责对一个恒星进行检索；

4)各线程检索出落入圆形视场内的恒星，在圆形视场内检索出落入CCD成像面的恒星。

2.如权利要求1所述的方法，其特征在于，所述步骤1)经纬度划分法是按照一定的经纬度间隔，用赤经圈和赤纬圈将天球分成不重叠的区域，每个区域内的恒星独立存储，并按照各个区域所在经度和纬度范围进行编号。

3.如权利要求2所述的方法，其特征在于，对同一纬度区间内的星区按照经度递增的顺序合并存储到一个星表文件，并且建立索引。

4.如权利要求3所述的方法，其特征在于，所述索引记录了各星区所在经度范围以及星区在星表文件中的恒星首地址、偏移长度。

5.如权利要求1所述的方法，其特征在于，所述步骤2)通过视场在天球的圆形投影面与赤经线构造球面三角形，通过球面几何关系计算出经纬度边界。

6.如权利要求4所述的方法，其特征在于，所述构造球面三角形方法为：以O为视轴与天球面的交点，过O点的经线与视场区域交与P₁、P₂两点，视场区域与天球两条相切经线切于Q₁、Q₂两点；所述球面几何关系为：

$\left\{\begin{array}{c}\cos (90-{\mathrm{\α}}_1)=\cos \mathrm{\θ}\cos (90-{\mathrm{\α}}_2)\\ \sin \mathrm{\θ}=\sin (90-{\mathrm{\α}}_1)\sin \mathrm{\Δ\δ}\end{array}\right.$

θ为视场半张角，α₁、α₂分别为O和Q₂的纬度，Δδ为Q₂和O的经度差。

7.如权利要求1所述的方法，其特征在于，所述步骤4)判断恒星落入圆形视场所在区域的方法是：计算出该恒星到观测点的连线与视轴的夹角，比较该夹角与视场张角的大小关系，如果该夹角大于视场张角，该恒星位于视场外，否则，恒星落入圆形视场所在区域。

8.如权利要求1所述的方法，其特征在于，所述步骤4)判断恒星坐标落入CCD像平面上的方法是：该恒星坐标在像平面的坐标范围内。

Images