地月圈层空间立体网格技术及其编码转换方法

张宗佩，万 刚，曹雪峰，李 锋，刘 婧

(信息工程大学地理空间信息学院，河南 郑州450052)

一、引 言

1959年前苏联发射的月球1号“梦想”成功进行近月飞行，标志着人类对月球探测活动的开始。自此之后，前苏联、美国竞相发射几十颗月球探测卫星，2004年我国也启动探月计划，目前已发射3颗对月探测卫星。人类的足迹早已涉足整个地月空间，地月空间中存在各种物质、现象，其对人类探月活动都有一定程度的影响。然而，目前缺乏一个统一的时空数据组织方法表达地月空间中存在的物质、现象及探月卫星运动过程。本文基于地球圈层空间立体网格理论［1］建立地球、月球圈层空间立体网格之间的联系，将地月立体空间统一到一个时空数据组织模型中，着重研究地月圈层空间立体网格基准面选择、编码转换及其涉及的地月空间坐标转换方法和时间基准，为在地月立体空间内统一描述各种实体属性、过程、现象建立严格的数学模型。

地月空间立体网格是地球圈层空间立体网格衍生出的新概念，是地球和月球圈层空间立体网格的结合体，整体结构如图1所示，其具有不同于地球、月球圈层空间立体网格的特性。地月圈层空间立体网格新特性有:

1)地球圈层空间立体网格覆盖月球立体空间边界。

2)月球圈层空间立体网格嵌套于地球圈层空间立体网格。

3)地球、月球圈层空间立体网格分别随地球、月球自转和公转。

图1 地月圈层空间立体网格整体结构

二、地月圈层空间立体网格时空基准

地球、月球都在不停地自转、公转，为了准确描述地月立体空间中各种时空现象，需要明确建立地月圈层空间立体网格及其时间和空间基准。

1．时空基准

时间和坐标是描述空间运动的基础要素，为了准确描述任何地月立体空间中物质、现象及运动，必须明确其运动的时间基准、空间坐标的基准。

(1)时间基准

地月立体空间内存在的活动对时间的精度要求很高，目前常用的时间系统［2-3］包括恒星时系统、太阳时系统、世界时、原子时、协调世界时和力学时。不同的时间系统对单位时间的定义不同，应用范围也各不相同。

文中为了统一描述时间轴上地球、月球自转和公转运动变化过程，以IAU建议的J2000．0(2 451 545．0)为初始时刻，即地球动力学时(TDB)2000年1月1日11:59:27．816，以儒略历元记录时间轴信息，基本时间单位为儒略日。

(2)空间基准

为了仿真地月立体空间中各种要素，必须建立准确的空间基准，即天球坐标系。地月立体空间中常用天球坐标系包括地心天体坐标系、地心地球天球坐标系、月心月球天球坐标系、月心月球赤道坐标系和月固坐标系［2］，坐标系定义见表1。

表1 5种坐标系定义［2］

地心天体坐标系和月固坐标系都是随着星体自转而旋转的坐标系，是建立地月圈层空间立体网格的空间坐标基准;而地心地球天球坐标系、月心月球天球和赤道坐标系是建立前两种坐标系间联系的中介。

2．地月空间坐标变换关系

为了在地月圈层空间立体网格内准确描述空间实体、空间现象，必须明确5种坐标系间变换关系。为了保证计算精度需求，采用IAU建议的旋转参数实现坐标系之间旋转变换［4］，采用VSOP87历表计算的行星位置实现坐标系之间平移变换，从而实现5种坐标系间坐标变换。IAU以国际天球参考架ICRF［3，5］作为基准，文献［4］中指出ICRF和J2000地心天球平赤道坐标系(即表1中地心地球天球坐标系)只存在小于1角秒的旋转偏差，因此文中采用后者作为IAU旋转的参考基准，如图2所示。记J2000．0时刻为t0，当前时刻为t。在IAU参考基准中行星北极点P(α0，δ0)，其中α0为赤经，δ0为赤纬，Q为行星赤道与参考基准赤道交点，其赤经为90°+α0，B为行星本初子午线与赤道交点，W表示QB之间旋转角度，行星赤道与参考基准赤道夹角为90°－δ0。

图2 行星旋转参数的参考框架示意图(IAU)

(1)地心天体与地心天球坐标变换

根据IAU提供的地球旋转参数α0、δ0、W，其参数实际数学含义如下

式中，d为相对于初始时刻t0的儒略日差值，即d=t－t0;T为相对于初始时刻t0的儒略世纪数，即T=d/36 525，单位为度。根据旋转参数可计算地心天体到地心天球坐标旋转矩阵ME为

通过旋转矩阵ME可实现地心天体坐标系到地心天球坐标系转换，同理可计算地心天球坐标系到地心天体坐标系的旋转矩阵M－1E。

(2)地心地球天球与月心月球天球坐标变换

从表1中可发现月心月球天球坐标系与地心地球天球坐标系两者坐标轴朝向定义一致，只有坐标原点不同。根据VSOP87历表可计算当前时刻t月心在地心天球坐标系中的位置为(Xt，Yt，Zt)，由于月心月球天球坐标系是由地心天球坐标系平移得到的，因此两者之间只存在一个平移向量

(3)月心天球与月固坐标变换

由于月心天球坐标系与地心天球坐标系仅坐标原点定义不同，且与月固坐标之间仅存在旋转变换，因此IAU提供的旋转模型可采用月心月球天球坐标系作为月球旋转模型的基准坐标系。

根据IAU提供的月球旋转参数α0、δ0、W，其数学含义如下

式中，E1=125°．045－0．052 992 1d，d和T含义同式(1)，完整公式和参数En(n=2，…，13)数学含义详见文献［5］。

因此，从月固坐标到月心月球天球坐标旋转矩阵ML为

同理，可计算得到从月心月球天球坐标到月固坐标的旋转矩阵M－1L。

三、地月圈层空间立体网格和编码转换技术

针对地月圈层空间立体网格新特性，不难发现现有的地球、月球圈层空间立体网格理论已难以满足其新需求。需要借助于建立的时空基准和现有的地球、月球圈层空间立体网格理论，满足地月圈层空间立体网格的新需求。

1．地月圈层空间立体网格基准圈层面

地月圈层空间立体网格是地球圈层空间立体网格［1］和月球圈层空间立体网格［6］的结合体。文献［1］给出地球的各个层次的基准圈层面(BSS)和基于基准圈层面的网格剖分方法。由于地月圈层空间立体网格的新特性要求和圈层空间立体网格可拓展特性，可拓展地球、月球的BSS，各层次BSS距离地心、月心的径向距离见表2。

表2 地球、月球基准圈层面

表2中地球半径Re=6 378．137 km，月球半径RL=1 734．0 km，—符号表示无数据。地球圈层空间外围边界达到72Re，即459 225．864 km，地月之间距离约为380 000 km，满足地月圈层空间立体网格新特性;月球圈层空间外围边界达到16RL，即27 744 km，可满足探月卫星等航天活动应用需求。

2．地月圈层空间立体网格编码转换模型

在地球或月球圈层立体空间内存在的实体、现象及其过程，若其仅在地球圈层立体空间内，则只需通过地球圈层空间立体网格编码进行描述表达;若其存在于地球和月球两个圈层立体空间内，则需要根据实际应用采用地球或月球圈层空间立体网格编码进行描述表达，避免在不同网格中进行操作。因此，必须对地月圈层空间立体网格编码转换进行研究。

地月圈层空间立体网格编码转换整体流程如图3所示。以地球圈层网格编码转换到月球圈层网格编码为例论述编码转换流程:

1)根据地球圈层网格编码、Hilbert解码算法［6］和地球地理坐标与网格编码对应关系［1］，计算圈层网格单元(记为圈体)左下角点P在当前时刻t的地心天体坐标(Xe，Ye，Ze)。

图3 地月圈层空间立体编码转换流程

2)通过第二章第2节给出的坐标旋转、平移变换，计算点P在月固坐标系中坐标(XL，YL，ZL)，即

3)根据月固坐标与月球网格编码对应关系、Hilbert编码算法［6］，可计算得出点P对应的月球圈层空间立体网格编码。

同理，可实现从月球圈层网格编码转换到地球圈层网格编码，从而建立地月圈层空间立体网格编码转换模型。

3．地月圈层空间立体网格编码转换方法

地月圈层空间立体网格编码转换最重要的是确定转换后的剖分层次，因此，在进行地球、月球圈层空间立体网格编码转换时，必须明确采用什么准则确定转换后编码的剖分层次。文中提出3种准则进行地月圈层空间立体网格编码转换，分别是等角度、等径向间距和等体积，并对以3种准则转换前后的圈体单元粒度和转换效率进行分析。

(1)等角度准则

等角度是指转换前圈体单元对应的球心角等于转换后圈体单元对应的球心角。根据地球、月球圈层空间立体网格剖分原理可知，保证相同的球心角等价于相同的剖分等级。因此，只要给出转换前圈体编码，就能确定转换后剖分等级，从而计算出转换后圈体编码。

由于地月之间距离大约为380 000 km，则仅在地球5-BSS和6-BSS之间的圈层体(最外层圈层体)才可能存在于月球圈层网格重叠的区域。月球圈层体最小径向间距为1000 km，地球最外层圈层体径向间距为40Re，转换前后径向间距比约为1∶255，则体积比约为1∶2553，空间粒度差距悬殊，但转换效率高。

(2)等径向间距准则

等径向间距是指转换前后的圈体单元在径向方向上间距相等，即

式中，ΔRL为月球圈层体径向间距;LevelL为月球剖分等级;ΔRE为地球圈层体径向间距;LevelE为地球剖分等级。将圈体单元近似抽象为一个长方体单元，根据体积解算公式，可得出转换前后地月圈体单元最小体积比约为1∶200。采用等径向间距准则进行编码转换，空间粒度差距较小，运算效率相比于等角度要低。

(3)等体积准则

等体积是指转换前后圈体单元的体积相等，即保持基本一致的空间粒度。为了便于计算，将圈格单元(如图4所示)简化为一个长方体，根据圈层网格剖分规则，得到近似等式为

式中，rE、rL分别为圈体单元左下角点与地心、月心的距离。由于sin 1/x与1/x正比的特性，将式(8)简化为

通过式(9)可利用近似等体积特性求解转换后圈层空间立体网格的近似剖分等级，从而实现地月圈层空间立体网格编码转换。此方法基本能实现空间粒度不变，但运算复杂，效率相较于前两者都低。

图4 圈格单元几何结构

四、试验与分析

根据地月圈层空间立体网格的时空基准和地月圈层空间立体网格剖分、编码算法，构建随时间绕太阳公转、自转的地月圈层空间立体网格，如图5所示。

图5 地月圈层空间立体网格

采用C语言实现3种转换准则下地月圈层空间立体网格编码转换算法，试验平台为Intel i5-4200M 2．5 GHz CPU、4 GB内存。计算10 000次从月球圈层空间立体网格编码到地球圈层空间立体网格编码消耗的时间，如图6所示。从图中可看出，3种转换准则效率从高到低依次为等角度、等径向间距、等体积。

图6 3种转换准则转换效率对比

五、结束语

本文详细论述了地月圈层空间立体网格所需时间基准和空间基准，并对地月空间坐标转换进行了深入研究，分析了地月圈层空间立体网格新特性，建立了地月圈层空间立体网格和编码转换模型，并分析比较了3种不同准则下编码转换效率和空间粒度变化情况。

［1］ 曹雪峰．地球圈层空间网格理论与算法研究［D］．郑州:信息工程大学，2012．

［2］ 马高峰．地月参考系及其转换研究［D］．郑州:信息工程大学，2005．

［3］ KOVALEVSKY J，SEIDELMANN P K．Fundamentals of Astrometry［M］．Cambrige:Cambrige University Press，2004．

［4］ ARCHINAL B A，A’HEARN M F，BOWELL E，et al．Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements:2009［J］．Celestial Mechanics and Dynamical Astronomy，2011(109):101-135．

［5］ PETIT G，LUZUM B．IERSTechnical Note No．36［R］．Frankfurt:IERSConventions，2010．

［6］ 李晨阳，张杨，冯玉才．N维Hilbert编码的计算［J］．计算机辅助设计与图形学学报，2006，18(7):1032-1038．

［7］ 韩阳，万刚，曹雪峰．混合式全球网格划分方法及编码研究［J］．测绘科学，2009，34(2):136-138．

［8］ 吴立新，余接情．基于球体退化八叉树的全球三维网格与变形特征［J］．地理与地理信息科学，2009，25(1):1-4．

［9］ MOKBEL F，AREFWG，KAMEL L．Analysis of Multidimensional Space-filling Curves［J］．GeoInformatica，2003，7(3):179-209．

［10］ 余接情，吴立新．球体退化八叉树网格编码与解码研究［J］．地理与地理信息科学，2009，25(1):5-9．

［11］ 蒋秉川，游熊，夏青．一种基于体素的三维地形可视化方法研究［J］．测绘通报，2013(3):46-49．