北斗静止轨道卫星信号盲区快速并行解算方法*

郭甲腾，吴立新，2，杨宜舟，李超岭，刘善军(.东北大学测绘遥感与数字矿山研究所，辽宁沈阳089；2.中国矿业大学物联网(感知矿山)研究中心，江苏徐州22008；.中国地质调查局发展研究中心，北京0007)

北斗静止轨道卫星信号盲区快速并行解算方法*

郭甲腾1，吴立新1，2，杨宜舟1，李超岭3，刘善军1
(1.东北大学测绘遥感与数字矿山研究所，辽宁沈阳110819；2.中国矿业大学物联网(感知矿山)研究中心，江苏徐州221008；3.中国地质调查局发展研究中心，北京100037)

北斗系统静止轨道卫星信号盲区解算方法复杂，串行计算耗费时间长，须在并行环境下利用更多的计算资源进行北斗盲区的快速解算。在分析北斗盲区解算原理与算法并行特征基础上，提出动态盲区影响域的并行解算方法。以栅格单元为并行粒度进行任务划分，实现了北斗盲区的高效并行解算。运用全国范围59景数字高程模型数据，利用8进程进行盲区并行解算，耗费时间约为5小时。实验测试结果表明:算法的并行效率随着进程数的增加有所衰减，但稳定在96%以上。其程序中间件已集成应用于高性能地理信息平台，应用效果良好。

北斗系统；并行计算；盲区解算；高性能地理信息系统；数字高程模型

北斗卫星导航系统(简称北斗系统)的空间段由5颗静止轨道卫星(北斗G4，G1，G3，G6与G5)和30颗非静止轨道卫星(中圆轨道卫星和倾斜地球同步轨道卫星)组成［1-3］。组成北斗系统的三种卫星的功能特点分别为:地球静止轨道卫星具备有源、无源、短报文通信三种服务功能［1-3］；中圆轨道卫星和倾斜地球同步轨道卫星只具有无源定位、导航和授时功能［1］。北斗系统中由于有了地球静止轨道卫星，才保证了有源定位、导航和授时服务以及短报文通信和位置报告功能。北斗静止轨道卫星的短报文通信功能可为采用北斗技术的户外移动设备、监测设备或野外工作人员提供短信播发服务，尤其针对手机信号盲区或无其他通信条件的地区(如无人区、高海拔地区等)，能在生命安全保障、生产物资调度、人员态势管理、日常业务沟通等方面发挥关键作用。目前，北斗卫星技术已正式投入民用，北斗卫星信号接收传感器也已日益普及［4］。由于北斗导航系统的优势，其被交通运输［5］、基础测绘［6］、工程勘测［7］等领域积极探索应用，尤其体现在应急救援［7-8］、物资保障［9］、信息推送［10-11］等方面。然而，北斗静止轨道卫星受地球表面弯曲、山体或障碍物遮挡等因素的影响，卫星终端设备在个别区域通信报文信号很弱(或称盲区)，制约着北斗通信的正常进行。因此，预先解算北斗系统信号盲区(或可视域)，能为基于北斗卫星应用技术的野外应用和管理工作提供重要的参考依据。由于北斗静止轨道卫星信号盲区解算(下文简称盲区解算)与地形直接相关，因此采用的基础数据是数字高程模型(Digital Elevation Model，DEM)数据。由于栅格DEM数据量庞大并且盲区解算复杂，采用串行算法解算效率较低。据前期经验，仅依据中国区域3″(90m)分辨率DEM的盲区解算就需花费数月，分辨率更高时则耗时更长，无法满足应急需求。近年计算机计算能力迅速提高，尤其是多核计算机和并行计算机开始普及，如何利用这些新的计算资源来解决更复杂的实际应用，已成为目前科学计算领域十分关注的问题［12］。在信息领域，并行计算为协同利用更多的计算资源提供了新的手段，实现了利用多个计算节点的计算资源来提高计算能力［13］。因此，郭甲腾等利用北斗静止轨道卫星数字化寻星原理［14］，引入地面点动态盲区影响域，提出了北斗信号盲区并行化解算方法，开发了并行计算中间件，将其集成到高性能地理计算平台(High Performance Geographic Information System，HiGIS)中，并基于SRTM 90m分辨率数据对全国范围北斗静止轨道卫星信号盲区进行了解算，大幅提高了计算效率。

1 北斗静止轨道卫星寻星原理

盲区解算中各参数如图1［14］所示:地球椭球长半轴a=6 378 137m；地球椭球短半轴b=6 356 752m；卫星s距地高度为H；卫星s的经度为ω；地面任意点p的海拔高度为h；p点经度为λ；p点纬度为φ；p点至卫星s的视域仰角为β；p点至卫星s的方位角为α；选定xoy平面为赤道面；x轴指向起始经度；z轴与旋转椭球体自转轴保持一致；e点为p点沿法线方向投影至椭球切平面处的交点；g点为过p点的法线向量与赤道面的交点；r为e点至原点o的距离；过p点的法截面与切平面的交线为ef；过p点且平行于切平面与os相交于q点；t为切平面与过p点的子午面交汇处。

地面任意点到卫星的距离为:

其中，gt为g点到t点的长度，et为e点到t点的长度，gs为g点到s点的长度，og为o点到g点的长度。

地面任意点到卫星的视域仰角为:

其中，θ为∠pgs，δ为∠psg。

图1 北斗卫星寻星参数几何关系图［14］Fig.1 Geometric relationship of parameters for finding BeiDou satellites

对于地面任意栅格单元p(下文简称目标栅格)，首先解算其与卫星位置连线的方位角，并计算出在方向上某一区域z(z为目标栅格的信号盲区影响域)内可能遮挡目标栅格卫星信号的栅格单元集S，并依次解算S中各栅格单元(下文简称遮挡栅格)针对目标栅格的信号遮挡状态。目标栅格卫星信号的遮挡栅格计算在区域z内均为独立计算，不同目标点的z值不同且互不影响，因此各栅格单元可并行地进行盲区计算。综上所述，实现盲区并行化解算的难点是目标栅格单元的盲区影响域z的计算。

2 运用DEM数据的北斗静止轨道卫星盲区并行解算

2.1 北斗信号盲区影响域分析

与北斗静止轨道卫星连线方位角的地面点呈线状分布(在DEM数据中为一组呈线性的栅格单元集合)，所以求解位置p的北斗信号状态(可视或盲区)需要方位角的其他所有栅格单元均参与计算。由于地球实际环境的影响，如地球曲率、地球高程差临界值ΔM(地球上海拔最高点与最低点最大相差约为20 000m，在中国区域内约为9200m)，所以目标栅格的盲区影响域应位于某一范围内，即不需要所有方位角的DEM栅格单元均参与栅格单元p的信号状态解算。p的卫星视域仰角与其所在的经纬度位置相关，p的经纬度值与卫星的经纬度值相差越大，则p的β越小［11］，盲区影响域范围越大。若d在某区域内的高程较大(如图2所示)，则其可在离p较远的位置对p的信号进行遮挡。若某区域内位置d高程最大，位置p高程最小，则两点的高程差Δm为该区域内的高程差临界值。因此，p的影响域范围大小与卫星的视域仰角及地球高程差临界值Δm相关，影响域半径为:

图2 北斗盲区影响域Fig.2 Influence scope of BeiDou signals blind area

信号盲区影响域的选择结果会直接影响到盲区计算的效率。由于中国大部分区域内局部范围的Δm远小于ΔM，若所有局部区域都采用中国区域的高差临界值ΔM来计算影响域范围z(ΔM)，将扩大局部区域的盲区影响域范围，从而使信号盲区的计算量急剧增加，特别是对处于北斗静止卫星可视边缘区域的影响尤为明显。因此，需在局部区域内计算位置p更为准确的盲区影响域的高程差临界值Δm(Δm＜ΔM)。所以针对中国范围，提出了局部区域内信号盲区影响域的计算方法。

首先采用中国区域内的高差极大值ΔM，根据式(4)获取以p为中心的z(ΔM)范围内数据(如图3所示)，遍历该区域内海拔最高栅格单元与海拔最低栅格单元的高差Δm，再以Δm为高差临界值，根据式(4)重新获取p的影响域z(Δm) (如图3所示)，通过z(Δm)逼近p实际的盲区影响域范围，以达到降低计算量的目的。

图3 p的盲区影响域Fig.3 BeiDou signals blind area of p

据上所述，被遮挡点p的信号影响域范围随着经纬度位置动态变化，其盲区计算量也不一样。而空间相邻的被遮挡点影响域大小基本相同，其盲区计算量差异较小。因此，在基于DEM的北斗盲区并行解算过程中，为确保各计算节点负载均衡，应将DEM数据中空间相邻的目标栅格分配至不同的计算节点进行解算。

2.2 信号盲区解算方法分析

运用上文所述盲区影响域计算方法，先解算DEM数据中所有栅格点的盲区影响域，在此基础上实现北斗盲区的解算(如图4所示)。p在局部范围内卫星信号被遮挡的区域(盲区影响域)大小为n×n(n=Re为栅格的分辨率)。

图4 北斗信号盲区解算Fig.4 Calculation of BeiDou signal blind area

1)计算p与北斗静止卫星连线的方位角α方向上且在z(Δm)内可能会对p信号遮挡的栅格单元集s(北斗静止轨道卫星在赤道上方，所以对于北半球，s在p的南方，而对于南半球，s在p的北方)。

2)p高程为h，相邻栅格点之间的距离为L，p的仰角为β，则距离p有m(1＜m＜n)个栅格单元位置的高程为hm，要遮挡p需满足的条件为:

若在大区域范围内则需要顾及地球曲率对高程的影响，表示为Δh=，其中，S为两点之间的距离，R为地球半径。

3)重复步骤2，依次计算s中各栅格单元对p的北斗信号遮挡状态，若发生遮挡，则结束计算，若不发生遮挡，则继续计算s中下一个目标栅格。

2.3 北斗信号盲区并行解算特点分析

基于DEM数据解算北斗信号盲区的过程，如图5所示，按照栅格A、B、C、D的顺序依次获取影响域并分别进行计算。北斗盲区计算属于计算密集型与数据密集型算法，将各栅格点与周围某一相关区域的点集进行计算，根据地理并行计算算法分类的方法［13］可知该算法属于区域计算，通过数据并行的方式即可实现盲区解算的并行化。虽然DEM中栅格单元的盲区影响域大小随着其经纬度的不同而动态变化，但空间相邻的栅格点的影响域范围大小基本相同，即相邻栅格单元中参与盲区计算的栅格点数目基本相等，所以基于DEM盲区计算具有以下特点:

1)将空间相邻的栅格单元分配至不同的计算节点，则各计算节点参与计算的栅格单元数量均衡。

2)由于空间相邻的栅格单元的盲区影响域大小基本相同(盲区影响域中参与计算的栅格数目基本相等)，将空间相邻的栅格单元分配至不同的计算节点，则各计算节点的计算量均衡。

图5 北斗信号盲区解算任务分布Fig.5 Task distribution for calculation of BeiDou signals blind area

现有并行编程模型可分为三类［15］:消息传递模型、共享存储模型和数据并行模型。共享存储模型属于线程级细粒度并行，数据共享存储、操作异步、数据隐式分配，程序可扩展性较差；数据并行模型属于进程级细粒度并行，数据共享存储、操作松散同步、数据半隐式分配，程序扩展性一般；而消息传递模型是目前使用最为广泛的并行模型，具有两大优势:

1)消息传递程序具有高度的可移植性与可扩展性，理论上可在任何并行机上执行，即不需要特殊的硬件支持。

2)允许用户显式地控制并行程序中每个进程内存的使用，为编程人员实现高性能计算提供便利。

消息传递模型的并行程序开发模式通常有两种:单程序多数据(Single Program Multiple Data，SPMD)模式和多程序多数据(Multiple Program Multiple Data，MPMD)模式。根据北斗盲区的解算特点，仅通过数据划分即可实现任务分解，或是通过任务量的划分也可实现任务的分解，且能保持北斗盲区并行计算过程中的负载均衡与任务均衡，即各进程对已分配任务执行相同的算法，对不同的数据集协同进行计算。因此，采用SPMD模式设计并行算法，实现北斗静止轨道卫星信号盲区的并行计算。

图6 北斗盲区并行解算流程Fig.6 Calculation process of BeiDou signals blind area

2.4 北斗信号盲区并行解算方法设计

根据上文分析可知，北斗盲区计算属于计算密集型算法(I/O与计算相比所耗费的时间较少)，算法以栅格点为单位进行解算，若直接通过数据划分(条带划分、块划分)，可能导致各节点的计算量不均衡(如6001×6001的DEM数据块分配给8个进程，若以条带划分会使其中7个进程计算量为750行，另外1个进程计算量为751行，使进程间差距6001个栅格点的计算量)，而通过以栅格点为粒度的任务划分，各进程间最多相差一个栅格点的计算量。由于解算过程中计算时间远大于DEM数据的I/O时间，因此采用任务划分的方法实现盲区计算的并行化，即各进程的DEM数据全冗余，通过任务量的划分方式实现并行任务的分解，如图6所示，将任务规划分至不同的计算节点实现北斗盲区解算的并行化。

针对北斗盲区计算，任务划分的粒度以栅格为单位，以保持各个进程的计算量均衡。任务划分的步骤如下:

1)待划分的DEM栅格数据集为Rs，其大小为M×N(DEM数据包含的栅格数目)，其中M为横向栅格数目，N为纵向栅格数目，DEM中每个栅格的位置标示为Lmn(lm，ln)；

2)参与栅格盲区计算的进程总数P={p1， p2，p3，……，pi}；

3)DEM数据是以行为单位组，以栅格单元为单位进行分配，若栅格单元的所在行号n与进程总数p取余得到pt，则以pt作为计算进程组的第一进程，对所有进程进行排序，排序后的结果为P= {pt，pt+1，…，pi，p1，…，pt-1}，如图7所示，通过此方法可将相邻的栅格单元分配至不同的计算节点；

4)计算进程组分配顺序后，以栅格单元的lm值为参考进行轮转分配；

5)重复步骤3、步骤4，直到所有的像素点分配完毕。

图7 任务划分Fig.7 Task partitioning

3 实验

3.1 测试环境

为测试方法的正确性与计算效率，搭建了一个小型集群的测试环境。该环境由2个高性能计算节点及1个存储节点(8T)构成，集群中各节点的CPU为Intel(R)Xeon(R)E5405(2.0GHz，64位)，操作系统为Redhat 6.6(64位)，消息传递采用OPENMPI并行库，节点间采用高性能交换机进行通信。

3.2 测试数据

实验采用的DEM数据为SRTM 90m分辨率(分辨率为3″)数据，每景数据的大小约为6001× 6001，数据量约为68MB。而中国区域内约包含59景数据，数据量约为4GB。5颗地球静止轨道卫星的参数见表1。

3.3 测试结果

为验证盲区解算方法的正确性，选择湖北、四川、陕西地区的交界处(精度范围为105°E～110°E，30°N～35°N的区域，该区域DEM数据大小为6001×6001)进行了信号盲区的并行解算。实验区域位于长江中游地区，以山区为主，有少量平原。经与中国地质调查局前期开展的信号盲区串行解算实验结果对比发现，实验区信号盲区基本一致，如图8所示。

表1 北斗同步轨道卫星信息Tab.1 Information of BeiDou geostationary orbit satellites

由北斗盲区并行解算性能(见表2)，可知:在单进程下解算上述区域的时间为1 092.69s，在8进程下并行解算的时间为142.65s，并行计算的最大加速比为7.66(并行计算比与进程数的关系如图9所示)，并行效率稳定在96%以上(并行效率与进程数的关系如图10所示)。利用8进程对全国范围内的59景DEM数据进行北斗卫星信号盲区并行解算，耗费时间约为5h。

图8 北斗盲区计算结果显示(a.中国北斗盲区矢量图；b.北斗盲区三维显示) Fig.8 Visualization of result of BeiDou signal blind area (a.Vector display of BeiDou signals blind area of China；b.Three-dimensional display of BeiDou signals blind area)

表2 北斗盲区并行解算性能Tab.2 Parallel computing performance ofBeiDou signal blind area

图9 北斗盲区并行解算加速比Fig.9 Speedup of BeiDou signal blind area calculation

同时，将北斗盲区并行计算的方法封装为中间件，并以工具集成应用到HiGIS平台中。在HiGIS平台中利用此方法进行计算，北斗盲区的解算结果与本文环境下的结果相同。

4 结论

针对北斗静止轨道卫星全国信号盲区并行解算的方法进行了探索与研究，引入了北斗信号盲区影响域的确定原理，提出了北斗信号盲区的快速求解方法，并在此基础上实现了盲区解算的并行化。由于测试实验环境为单机多核，当进程数增加时，并行加速比会随着进程数的增加而有所衰减，但基本稳定在96%以上，并行加速比达到7.6以上。若在网络通信性能较弱的集群上进行测试，则需考虑I/O的影响。实验结果表明，北斗静止轨道卫星盲区并行解算方法可行。同时，基于本文方法实现的程序中间件已集成部署在HiGIS平台中，有望应用于高性能计算环境下基于更高精度DEM数据的北斗静止轨道卫星信号盲区快速解算。

图10 北斗盲区解算并行效率Fig.10 Parallel efficiency of BeiDou signal blind area calculation

References)

［1］BDS-OS-PS-1.0.北斗卫星导航系统公开服务性能规范［S］.北京:中国卫星导航系统管理办公室，2013.BDS-OS-PS-1.0.BeiDou Navigation satellite system open service performance standard［S］.Beijing:China Satellite Navigation Office，2013.(in Chinese)

［2］Yang Y X，Li JL，Wang A B，et al.Preliminary assessment of the navigation and positioning performance of BeiDou regional navigation satellite system［J］.Science China:Earth Sciences，2014，57(1):144-152.

［3］杨元喜，李金龙，徐君毅，等.中国北斗卫星导航系统对全球PNT用户的贡献［J］.科学通报，2011，56(21): 1734-1740.YANG Yuanxi，Li Jinlong，Xu Junyi，et al.Contribution of the compass satellite navigation system to global PNT users［J］.Chinese Science Bulletin，2011，56(21):1734-1740.(in Chinese)

［4］唐金元.北斗卫星导航区域系统发展应用综述［J］.全球定位系统，2013，38(5):47-52.TANG Jinyuan.Developing and applying analysis of BeiDou navigation satellites regional system［J］.GNSS World of China，2013，38(5):47-52.(in Chinese)

［5］万琪俊.基于北斗导航卫星移动定位系统研究与实现［D］.广州:华南理工大学，2013.WAN Qijun.Research and implementation ofmobile location system of BeiDou navigation satellite base on android［D］.Guangzhou:South China University of Technology，2013.(in Chinese)

［6］王家耀.测绘导航与地理信息科学技术的进展［J］.测绘科学技术学报，2014，31(5):441-449.WANG Jiayao.The progress of geomatics navigation and geographic information science and technology［J］.Journal of Geomatics Science and Technology，2014，31(5):441-449.(in Chinese)

［7］温静，汪大明，孟月玥，等.北斗卫星导航系统在地质调查领域应用综述［J］.地质力学学报，2012，18(3): 213-223.WEN Jing，WANG Daming，MENG Yueyue，et al.Application of BeiDou navigation satellite system to geological survey［J］.Journal of Geomechanics，2012，18(3):213-223.(in Chinese)

［8］石小亚，王占昌，伍锦程.北斗卫星技术在青海玉树地质灾害详查中的示范应用［J］.地质力学学报，2012，18(3):277-281.SHIXiaoya，WANG Zhanchang，WU Jincheng.Demonstration application of BeiDou satellite technology in detailed geological disaster investigation in Yushu，Qinghai province［J］.Journal of Geomechanics，2012，18(3):277-281.(in Chinese)

［9］李小杰，叶成明，郑宝锋.基于北斗的地质调查安全生产保障系统设计与实现［J］.安全与环境工程，2015，22(3):136-139.LI Xiaojie，YE Chengming，ZHEN Baofeng.Design and implementation of the guarantee system for geological survey safety production based on BeiDou satellite navigation［J］.Safety and Environmental Engineering，2015，22(3):136-139.(in Chinese)

［10］何凯涛，李志忠，汪大明.基于遥感卫星和北斗卫星的野外地质调查服务与管理系统设计综述［J］.地质力学学报，2012，18(3):203-212. HE Kaitao，LI Zhizhong，WANG Daming.Overview on the design of the service and management system for field geological survey based on the remote sensing and BeiDou satellites［J］.Journal of Geomechanics，2012，18(3): 203-212.(in Chinese)

［11］伍锦程，王占昌，石小亚.西北区野外地质调查中北斗卫星通讯与定位技术应用研究［J］.地质力学学报，2012，18(3):282-287.WU Jincheng，WANG Zhanchang，SHIXiaoya.Application of BeiDou satellite communication and positioning technology to the field geological survey in northwestern china［J］.Journal of Geomechanics，2012，18(3):282-287.(in Chinese)

［12］赵春宇.高性能并行GIS中矢量空间数据存取与处理关键技术研究［D］.武汉:武汉大学，2006.ZHAO chunyu.Studying on the technologies of storage and processing of spatial vector data in high-performance parallel GIS［D］.Wuhan:Wuhan University，2006.(in Chinese)

［13］吴立新，杨宜舟，秦承志，等.面向新型硬件构架的新一代GIS基础并行算法研究［J］.地理与地理信息科学，2013，29(4):1-8.WU Lixin，YANG Yizhou，QIN Chengzhi，et al.On basic geographic algorithms of new generation GIS for new hardware architectures［J］.Geography and Geo-Information Science，2013，29(4):1-8.(in Chinese)

［14］王占昌，王飞，伍锦程，等.一种静止轨道地球同步卫星数字寻星方法:中国，CN201410020373［P］.2014-01-06.WANG Zhanchang，WANG Fei，WU Jincheng，et al.A method of digital finding satellite for geosynchronous satellites:China，CN201410020373［P］.2014-01-06.(in Chinese)

［15］Xavier C，Iyengar SS.Introduction to parallel algorithms［M］.USA:John Wiley＆Sons，Inc.，2004.

Parallel signal blind area computing of BeiDou geostationary orbit satellite

GUO Jiateng1，WU Lixin1，2，YANG Yizhou1，LIChaoling3，LIU Shanjun1
(1.Institute of Geo-informatics＆Digital Mine，Northeastern University，Shenyang 110819，China；2.IoT/Perception Mine Research Center，China University of Mining＆Technology，Xuzhou 221008，China；3.Development and Research Center，China Geological Survey，Beijing 100037，China)

In order to improve the computing efficiency of signal blind area for BeiDou geostationary orbit satellites，it requiresmore computing resources in parallel computing environment.According to the analysis on parallel characteristics of the existing resolving principle and algorithm，a novelalgorithm was designed for computing signal blind area of BeiDou geostationary orbitsatellitesbased on elevation angle related dynamic blind influence scope.Based on this algorithm，the parallel computing tasks of signal blind area was divided and arranged by raster data partition.An experimentwas taken for computing the blind area in china based on 59 scene digital elevation model data and 8 processes，and spent about 5 hours.The result shows that the parallel efficiency of the algorithm decreaseswith the increasing process count，but stablymaintains atmore than 96%.Themethod has been integrated as a tool in high performance geographic information system and can obtain good application effect.

BeiDou system；parallel computing；blind area calculation；high performance geographic information system；digital elevation model

TP391

A

1001-2486(2015)05-047-07

10.11887/j.cn.201505008

http://journal.nudt.edu.cn

2015-06-29

国家863计划资助项目(2011AA120302)；国家自然科学基金资助项目(41001228)；中央高校基本科研业务费资助项目(N140104002)；辽宁省自然科学基金资助项目(2015020581)

郭甲腾(1980—)，男，安徽桐城人，讲师，博士，E-mail:guojiateng@mail.neu.edu.cn