随着传感器及探测技术的迅猛发展,地球空间信息的获取、处理和应用已步入大数据时代^[1]。对地观测数据海量、广域、多源、异构、多尺度等特点使得地球大数据的表征、汇聚和处理愈加困难。全球离散格网系统（discrete global grid systems,DGGS）是对整个地球空间进行递归剖分形成的多分辨率离散地球参考模型^[2],有助于建立以空间位置为主键的数据关联新模式,有望成为支持地球大数据融合分析的优选方案^[3-4]。

基于多面体投影是构建DGGS的有效方法。直接在球面进行规则格网剖分较为困难,采用特定方法递归剖分理想多面体表面,再映射至地球表面,即可得到较为规则的球面格网。多面体投影能在满足等积或等角性质的情况下,另一性质的变形较小,使构建格网的几何性质较为优异。此外,多面体格网与球面格网具有相同的拓扑关系,原本在球面的部分运算操作可转换至多面体表面,有效提高了计算效率。

DGGS的构建主要涉及基础多面体、投影方法、多面体定位、剖分方案与编码5个设计选项^[5],其中基础多面体与投影方法是影响格网系统单元形变的主要因素。在基础多面体方面,由于柏拉图多面体规则的几何结构,现有绝大部分DGGS采用该类多面体来拟合地球^[6-9]。文献[10]系统性地统计分析了不同柏拉图多面体与投影方法对格网变形性质的影响,发现面数最多的正二十面体带来的畸变最小。目前正二十面体已成为构建DGGS的主要方案^[8,11-13],但其仍然存在与球面差异过大、格网单元分布不均匀的问题。为突破柏拉图多面体的束缚,文献[14]尝试采用卡特兰多面体中的120面体构建DGGS,相比正二十面体,该方案的角度变形大大减小,但其单个面均为不规则三角形,仅适于构造三角形格网系统,且计算较复杂。

在投影方面,单元等积性被认为是较为理想的格网特性,它能有效提高数据表达精度并减少数据冗余^[15],适于统计分析计算。多面体等积投影相关研究较少,现有柏拉图多面体DGGS主要采用Snyder等积多面体投影^[16],该投影解决了相邻投影面公共边经纬线断裂的问题,然而在投影边界处的角度变形较大,且逆投影需迭代求解,有损精度并影响效率。文献[17]提出了slice and dice多面体等积投影,简称Leeuwen投影,相比Snyder投影,角度变形更均匀,但仅用于柏拉图多面体,且未给出多面体到球面的逆投影公式。DGGS格网生成、数据分析与可视化等方面需频繁使用逆投影,研究表明其运算效率已成为限制DGGS应用效率的瓶颈之一^[18],因此有必要研究更高效的多面体逆投影运算。

综上,选择合适的多面体并设计高效的投影方法可以构建更加规则、均匀,且性能优异的DGGS。本文选择菱形三十面体作为基础多面体,虽然面数少于120面体,但其结构在DGGS应用中存在以下优势：（1）面数多于正二十面体,与地球拟合程度更高,且相邻两个面沿公共边镜像对称符合多面体投影基本要求；（2）每个面为大小相同的菱形,适合各种规则单元的剖分。尽管八面体或正二十面体可利用两个三角面组成菱形作为剖分及编码的基本结构^[7-8],但其仅是一种逻辑菱形结构,两个三角面不位于同一平面,增加了单元位置描述及格网相关运算的复杂度；（3）菱形面上的一组邻边是天然的二维坐标系,易于描述单元位置。结合上述结构特点,依据Leeuwen投影原理,本文推导了菱形三十面体正逆投影的解析表达式,并通过对比实验验证了该方案在格网变形与计算效率上的优势。

1 基本原理

Leeuwen投影是实现球面三角形到多面体三角面等积映射的一类方法^[17],其基本原理是：在某种分割策略下,球面三角形上任意无限小的单元均以相等的面积投影至平面上一无限小单元,该投影即满足等积性质。根据分割策略的不同,Leeuwen等积投影有不同实现方法,本文采用的是指向多面体顶点的大圆弧等积投影,该投影的分割策略为：从顶点出发分割球面三角形的大圆弧,投影到平面三角形为过对应顶点的一条直线,其中顶点指的是多面体及其球面多面体的顶点,如图1所示。基于上述分割策略,根据Leeuwen投影原理,可推导出该投影的等积条件,推导过程见文献[17]。

图  1  面向顶点的大圆弧投影示意图

Figure  1.  Diagram of Vertex-Oriented Great Circle Projection

下载: 全尺寸图片 幻灯片

如图1（a）所示,$P$为球面三角形$T$上任意一点,$l$为穿过顶点$B$且经过点$P$的大圆弧,$l$与球面三角形的边$AC$交于点$D$并将$T$分为区域$U$与$T-U$,点$P$将大圆弧$BD$分割为两部分,令弧长$\stackrel{̑}{BP}=x$,$\stackrel{̑}{PD}=y$；如图1（b）所示,$T^{\text{'}}$为多面体表面上的投影三角面,$P^{\text{'}}$为$P$在$T^{\text{'}}$上投影点,$l^{\text{'}}$为过平面三角形顶点$B^{\text{'}}$与$P^{\text{'}}$的一条直线,$l^{\text{'}}$与边$A^{\text{'}}{C}^{\text{'}}$交于点$D^{\text{'}}$并将$T^{\text{'}}$分割为区域$U^{\text{'}}$与$T^{\text{'}}-U^{\text{'}}$,线段$C^{\text{'}}{D}^{\text{'}}=u^{\text{'}}$,$A^{\text{'}}{D}^{\text{'}}=v^{\text{'}}$,点$P^{\text{'}}$将线段$B^{\text{'}}{D}^{\text{'}}$分割为$B^{\text{'}}{P}^{\text{'}}$和$P^{\text{'}}{D}^{\text{'}}$,定义线段长度$B^{\text{'}}{P}^{\text{'}}=x^{\text{'}}$,$P^{\text{'}}{D}^{\text{'}}=y^{\text{'}}$。指向多面体顶点的大圆弧等积投影条件如下：

1）球面三角形$T$的面积$A_T$与平面三角形$T^{\text{'}}$的面积$A_{T^{\text{'}}}$相等,即

2）球面三角形上区域$U$与$T$的面积比等于平面三角形上区域$U^{\text{'}}$与$T^{\text{'}}$的面积比,即

3）在大圆弧$l$与其投影直线$l^{\text{'}}$上,点$P$投影点$P^{\text{'}}$的位置满足：

上述3个条件中,式（1）和式（2）适用于所有基于Leeuwen投影原理的等积投影,式（3）是根据具体的分割策略推导等积微分条件而得。

2 菱形三十面体等积投影

菱形三十面体由30个全等菱形组成,共有32个顶点、60条边,每个面均为黄金菱形,即两条对角线的长度比为黄金比例（图2（a））。Leeuwen投影实现了球面三角形到平面三角形的等积映射,为将该方法推广至菱形三十面体,本文将三十面体的各个菱形面沿短对角线分为两个相等的等腰三角形作为菱形三十面体的基本投影面。考虑到菱形三十面体的对称特性,下文仅针对一个等腰三角形面进行投影算法的推导。

图  2  菱形三十面体正逆投影示意图

Figure  2.  Diagram of the Forward and Inverse Projections Between Rhombic Triacontahedron and Sphere

下载: 全尺寸图片 幻灯片

2.1   投影基本参数

多面体投影需先确定球面三角形与平面三角形的顶点坐标、边长等作为投影运算的基本参数。为计算菱形三十面体基本投影面,即等腰三角形的有关参数,将等积条件式（1）推广到整个多面体,等效于菱形三十面体表面积等于球面面积,由此可计算与单位球等表面积的菱形三十面体菱形面棱长$r$、长对角线$p$、短对角线$q$分别为：

为方便平面坐标运算,在图1（b）所示平面三角形$T^{\text{'}}$上构建二维平面直角坐标系：以$A^{\text{'}}{C}^{\text{'}}$为$x$轴,$A^{\text{'}}{C}^{\text{'}}$中点为原点,$B^{\text{'}}$位于$y$轴正方向上,易得3个顶点坐标分别为$A^{\text{'}}\left(-q/\mathrm{2,0}\right)$、$C^{\text{'}}\left(q/\mathrm{2,0}\right)$、$B^{\text{'}}\left(\mathrm{0},p/\mathrm{2}\right)$。

为计算球面菱形三十面体（图2（b））上顶点坐标、边的弧长与内角等参数,需先定义菱形三十面体与地球球体的位置关系。球面菱形三十面体共32个顶点,以数字0~31对其索引,在展开面上的具体标识如图3所示。

图  3  菱形三十面体顶点索引

Figure  3.  Vertex Indices of the Rhombic Triacontahedron

下载: 全尺寸图片 幻灯片

本文采取一种简单、易理解的定位方法：将0号顶点放置在北极,31号顶点放置于南极,1号顶点经度为36°E,可计算得菱形三十面体32个顶点的地理坐标,1~5、6~10、11~15顶点纬度分别为52.622 631 859°N、26.565 051 177°N、10.812 316 964°N,同纬度相邻顶点经度间隔72°,根据对称性,可得其他顶点地理坐标。该定位方法使得投影后经度为东经（或西经）0°、36°、72°、108°、154°、180°的10条经线均与菱形面的对角线或棱长重合,且关于赤道对称。在格网构建时,能使离散化的球面格网单元与地理特性建立关联,利于地球科学相关研究。在实际情况中也可根据应用需求,选择其他定位参数,并计算相应的顶点地理坐标。

根据顶点坐标,利用球面几何原理可计算得到球面三角形T的相关参数：$\mathrm{\angle }A=\mathrm{\angle }C=\gamma =\mathrm{\pi }/\mathrm{3}$,$\mathrm{\angle }B=\beta =\mathrm{2}\mathrm{\pi }/\mathrm{5}$,$\stackrel{̑}{BA}=\stackrel{̑}{BC}=\mathrm{0.652}\mathrm{358}\mathrm{139}$,$\stackrel{̑}{AC}=\mathrm{0.729}\mathrm{727}\mathrm{656}$。

2.2   正投影公式

给定球面三角形上一点地理坐标,解算其在平面三角形上投影坐标的过程为正投影或投影正算。根据等积投影条件,推导菱形三十面体等积投影正算的解析表达式。

1）根据等积条件式（2）确定投影点$P^{\text{'}}$所在投影直线$l^{\text{'}}$在平面三角形$T^{\text{'}}$上的位置。

依据单位球上球面三角形的面积公式可得$A_T=\beta +\mathrm{2}\gamma -\mathrm{\pi }=\mathrm{\pi }/\mathrm{15}$；设大圆弧$l$与$AB$、$AC$的夹角分别为$\rho$和$\delta$,则可得球面三角形$ABD$对应区域$U$的面积$A_U=\frac{\mathrm{\pi }}{\mathrm{15}}-(\rho +\delta +\frac{\mathrm{\pi }}{\mathrm{3}}-\mathrm{\pi })=\frac{\mathrm{11}\mathrm{\pi }}{\mathrm{15}}-\rho -\delta$。因此,可以推导在菱形三十面体等积投影中投影条件式（2）$\frac{A_U}{A_T}$的具体表达式为：

式中,由于点$P$的经纬度已知,可用叉乘公式解算$\rho$,然后在球面三角形$ABD$内根据球面余弦公式计算$\delta$,计算式为：

将$\rho$与$\delta$代入式（5）可得$\frac{u^{\text{'}}}{u^{\text{'}}+v^{\text{'}}}$,易计算直线$l^{\text{'}}$与$A^{\text{'}}{C}^{\text{'}}$的交点$D^{\text{'}}$坐标$\left(x_{D^{\text{'}}},y_{D^{\text{'}}}\right)$：

2）根据等积条件式（3）确定投影直线$l^{\text{'}}$上投影点$P^{\text{'}}$的位置。

已知点$B$、$P$的地理坐标,计算两点间弧长$x$；在球面三角形$ABD$内,借助余弦公式计算$\mathrm{c}\mathrm{o}\mathrm{s}(x+y)$：

根据等积条件式（3）,可得$\frac{x^{\text{'}}}{x^{\text{'}}+y^{\text{'}}}$,然后计算投影点$P^{\text{'}}$坐标$\left(x_{P^{\text{'}}},y_{P^{\text{'}}}\right)$：

2.3   逆投影公式

给定平面三角形上一点的平面坐标,解算其在球面三角形上投影点地理坐标的过程为逆投影或投影反算。同样根据等积投影条件,推导菱形三十面体等积投影反算的解析表达式。

1）求解球面上过投影点的大圆弧与球面三角形交点$D$的坐标。

对于逆投影,点$P^{\text{'}}$的平面坐标已知,易计算等积投影条件式（2）的比例系数$\frac{u^{\text{'}}}{u^{\text{'}}+v^{\text{'}}}$,然后根据式（5）计算$\rho +\delta$。

为解算$\rho$与$\delta$,根据球面三角形ABD内余弦公式$\mathrm{c}\mathrm{o}\mathrm{s}\delta =-\mathrm{c}\mathrm{o}\mathrm{s}\rho \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\rho \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\stackrel{̑}{AB}$,令$\tau =\rho +\delta$,并用$\tau -\rho$替换$\delta$,展开得：

两边除以$\mathrm{c}\mathrm{o}\mathrm{s}\rho$得：

即可解算出球面三角形$ABD$内角$\rho$与$\delta$,再根据余弦公式可求出$\mathrm{s}\mathrm{i}\mathrm{n}\stackrel{̑}{AD}$。

由于$\mathrm{s}\mathrm{i}\mathrm{n}\stackrel{̑}{AC}$已知,根据四元数的球面线性插值（spherical linear interpolation,SLERP）,得到点$D$的坐标计算式为：

式中,$q_D$、$q_C$、$q_A$分别为球面三角形对应点的四元数坐标。由于投影坐标求解无需旋转参数,因此只用点$A$、$C$的三维坐标即可计算出点$D$的三维坐标。

2）求解大圆弧$\stackrel{̑}{BD}$上投影点$P$的坐标。再次利用SLERP插值方法,大圆弧$\stackrel{̑}{BD}$上有：

如图1所示,$\stackrel{̑}{BD}=x+y$,$\stackrel{̑}{BP}=x$,$\stackrel{̑}{PD}=y$。为求解点$P$坐标,需计算上述弧长。根据投影基本参数,易得点$B$、$D$三维坐标,然后计算$x+y$；对于$x$,已知平面三角形内点坐标,计算$x^{\text{'}}$与$y^{\text{'}}$,然后根据式（3）,即可解算$x$。

将$x+y$、$x$与$y$代入式（13）计算$q_P$,即为投影点$P$坐标。可发现,上述运算过程均为解析表达式,无需Snyder逆投影的迭代过程。

3 实验与分析

为了验证本文所提方案的优势,设计以下实验：（1）投影角度变形对比实验,比较Leeuwen投影与Snyder投影应用于菱形三十面体与正二十面体的角度畸变；（2）格网生成实验,分别生成三角形、四边形和六边形三类规则单元格网系统；（3）格网生成效率对比实验,比较本文推导的逆投影与Snyder等积多面体逆投影对格网生成效率的影响。

3.1   等积投影角度变形对比实验

由于地图投影无法同时满足等积和等角性质,因此等积投影必然存在角度变形。为证明本文方案的优势,比较其与目前使用广泛且面数最多的正二十面体Snyder等积投影的角度变形。实验采用Tissot椭圆^[19]作为角度变形程度的通用衡量指标来计算角度变形值,其原理为：球面上任意无限小的圆投影至多面体表面为一椭圆。定义投影椭圆的长轴和短轴分别为$a$和$b$,角度变形值$d$的计算式为：

具体实验方法为：在球面三角形上,选择均匀分布的若干小圆（图4）,每个小圆上等间隔取若干离散点,投影至平面三角形,计算小圆上投影点与小圆中心投影点的距离,取距离中的最大值与最小值分别作为此小圆投影椭圆的长轴和短轴,然后利用式（14）解算小圆的角度变形$d$,最后统计所有小圆角度变形均值与标准差。

图  4  Tissot椭圆球面小圆采样示意图

Figure  4.  Sampled Small Circles of Tissot Ellipse on Sphere

下载: 全尺寸图片 幻灯片

经过实验统计,各投影角度变形的均值与标准差如表1所示。为更直观展示角度变形程度,将角度变形值可视化,按照1°间隔,将0°~20°的角度变形值渐变映射为白色至品红色,其中20°是由实测角度变形最大弧度值——正二十面体Snyder投影的0.31确定。可视化结果如图5所示,颜色越浅代表角度变形的程度越小,颜色越深代表角度变形的程度越大。

表  1  等积投影角度变形均值与标准差/rad

Table  1.  Average Value and Standard Deviation of the Angular Distortion for Different Equal-Area Projections/rad

投影类型	均值	标准差
正二十面体Snyder	0.166	0.055
菱形三十面体Leeuwen	0.082	0.023

下载: 导出CSV 

| 显示表格

图  5  等积投影角度变形可视化结果

Figure  5.  Visualization of the Angular Distortion for Different Equal-Area Projections

下载: 全尺寸图片 幻灯片

实验结果证明了本文设计方案在基础多面体与投影方法上的优势,分析实验结果可发现：（1）菱形三十面体Leeuwen投影的整体角度变形与二十面体Snyder投影相比明显减小,且分布更加均匀；相比二十面体Snyder投影,菱形三十面体Leeuwen投影的平均角度变形减小约51%,标准差减小约58%,原因主要在于菱形三十面体的面数更多,与球体的近似程度更高,因此变形更小。（2）由于投影策略不同,两种投影的变形分布也不同：Snyder投影是以三角面中心作为投影点,角度变形从三角面中心向外辐射,可以直观地看出,三角面内部3条投影边界处变形非常大,角度变形接近20°；菱形三十面体Leeuwen投影以菱形面的顶点为投影中心,每个面上角度变形从投影中心向外呈辐射状散发,离投影中心越远,变形越小,面中心处几乎无角度变形。

3.2   菱形三十面体格网生成实验

菱形三十面体的30个面为全等的菱形,其优势在于菱形面适合多种规则图形的剖分,包括三角形、菱形和六边形。本文基于菱形三十面体等积投影方法构建全球等积格网的基本流程如图6所示。

图  6  菱形三十面体等积全球离散格网系统构建流程图

Figure  6.  Flowchart of Constructing Rhombic Triacontahedron Equal-Area Grid Systems

下载: 全尺寸图片 幻灯片

确定多面体及与地球定位关系后,格网系统构建的下一步是确定单元类型及初始剖分单元。根据菱形三十面体结构,本文定义如下三类格网系统的初始单元,即第0层单元：30个菱形面为菱形格网系统第0层单元；每个菱形面沿对角线分为两个三角形,共60个三角面作为三角形格网系统第0层单元；由60个三角面组成的12个五边形作为六边形格网系统的第0层单元,五边形分别以图3中索引为0、6~10、21~25与31的12个顶点为中心。

较特别的是六边形格网,以格网剖分孔径4为例,即上下两层单元面积之比为4,六边形格网第1层单元分布如图7所示,共30个六边形和12个五边形。需要注意的是,六边形格网系统每层都是由六边形和12个五边形组成,五边形的面积为六边形的5/6。

图  7  孔径为4菱形三十面体六边形格网系统第1层单元

Figure  7.  The First Level of Aperture-4 Hexagonal Grid System Based on Rhombic Triacontahedron

下载: 全尺寸图片 幻灯片

图8展示了算法生成的菱形三十面体三角形、菱形、六边形球面等积格网的第三层,验证了本文所提等积投影算法能够很好地适用于三类规则单元格网系统构建。

图  8  三类规则单元的菱形三十面体等积格网第3层示例

Figure  8.  The Third Level of Three Different Types of Equal-Area Grids Based on Rhombic Triacontahedron

下载: 全尺寸图片 幻灯片

3.3   格网生成效率对比实验

在格网生成过程中,逆投影运算需频繁使用三角函数运算,是制约格网生成效率的主要因素。为进一步验证本文提出逆投影算法的效率优势,采用Visual C++ 2012开发工具分别实现了基于二十面体Snyder逆投影与本文提出菱形三十面体等积逆投影的格网生成算法。

为获取均匀分布的待投影单元集合,分别构建二十面体和菱形三十面体三角形格网,相邻两层单元的面积比均为4,计算三角形单元的顶点坐标作为待投影点,并投影至球面计算地理坐标。全部程序编译为Release版本,并在一台兼容机上进行测试,操作系统为Windows 7 x64 旗舰版+SP1,硬件配置为Intel Core i5-6500 CPU@3.20 GHz,8 GB RAM,KingSton 120 GB SSD。统计各算法在格网7~10层的运算时间,运算效率为单位时间内生成的三角形单元个数,每层测试10次取平均值,结果见表2。

表  2  格网生成实验效率统计结果

Table  2.  Efficiency of the Global Grid Generation for Icosahedron Snyder Grids and the Proposed Grids

层次	本文方案			二十面体Snyder			效率比
	单元数量	运行时间/ms	效率/m·s-1	单元数量	运行时间/ms	效率/m·s-1
7	983 040	760	1 293.47	327 680	997	328.67	3.94
8	3 932 160	3 003	1 309.41	1 310 720	4 021	325.97	4.02
9	15 728 640	12 047	1 305.61	5 242 880	15 901	329.72	3.96
10	62 914 560	47 840	1 315.10	20 971 520	6 2718	334.38	3.93

下载: 导出CSV 

| 显示表格

分析实验结果可知：（1）本文算法生成菱形三十面体等积格网的效率约为二十面体Snyder等积投影生成格网效率的4倍。除逆投影算法外,两者生成过程中的其他运算基本一致,均为整数或简单的浮点数运算,可见本文逆算法对格网生成效率的提升。主要原因在于本文算法求解过程均为解析表达,而Snyder逆投影涉及迭代运算,且需良好的初值,大大影响了运算效率。（2）影响格网生成效率的主要因素是逆投影算法的三角函数运算,由于两种算法均涉及难以避免的球面三角形参数解算,因此效率提升存在瓶颈。

4 结语

本文引入菱形三十面体提出了一种新型等积全球离散格网系统构建方法。根据等积投影条件,推导了菱形三十面体正逆等积投影的解析表达式,并利用逆投影算法构建了三类不同形状的等积全球格网。与现有方案相比,该方案构建的DGGS单元角度变形更小,有利于提高后续数据处理与表达的精度,且解决了现有多面体等积逆投影迭代导致格网构建效率低的问题。本文方案摆脱了当前利用柏拉图多面体构建DGGS的思维定势,为全球离散格网系统的发展及相关应用提供了新的思路。下一步需基于该方案继续研究数据建模、空间分析算法及DGGS与高性能计算的结合,为地球大数据的组织、处理与分析提供更优的解决方案。