球面四元三角网(quarternary triangular mesh，QTM)结构^[1]具有层次组织、全球连续和格网单元形状统一等特性，是目前较为基础、最有影响、应用最为普遍的球面离散格网剖分模型^[2]，在全球数据组织与管理^[3-7]、全球数据质量与地图综合^{[1, 8]}、全球地形可视化^[9-13]、全球多分辨率水淹模拟^[14-15]、全球影像数据检索^[16]、全球导航^[17]等领域均有广泛的应用。作为一种球面栅格数据结构，QTM对地理空间数据和信息的表达是建立在矢量数据向QTM地址码转换的基础之上的。对于不断更新的全球海量数据，地址码与经纬度的转换精度及效率已经成为影响系统应用的主要因素^[18]。

针对该问题，许多学者进行了深入研究。例如，Goodchild等^[4]提出了等三角投影法(equal-triangles projection, ETP)，利用ETP投影和距离比较的思想准确地对地址码和经纬度进行了转换，但其计算过程复杂，需要求幂和无理算术等操作，效率较低。Dutton^[19]提出了天顶正交(zenithal ortho triangular, ZOT)投影法，在ZOT空间利用Manhattan距离进行转换，将转换效率提高了一个量级，但该算法转换得到的地址码丧失了方向性，不利于计算机的邻近搜索。为解决转换速度和编码方向性之间的矛盾，部分学者提出了行列逼近法^[18]及三向互化算法^[20]，这两种算法的转换效率均与ZOT投影法处于同一量级，且生成的地址码具有方向性，但也产生了新的问题。行列逼近法在经纬度向地址码转换的过程中多出了半个格网的误差，三向互化算法则在格元边界处存在“交叉地带”。这种误差的存在为球面QTM模型空间应用分析与操作埋下了隐患。

为建立一种既保证转换效率与方向性（以适应于当今海量数据应用需求），又保证较高精度的QTM地址码与经纬度转换算法，本文在行列逼近法和ETP投影法的基础上进行了一系列改进，即建立行列号(i，j)与球面经纬度坐标(λ，φ)以及球面三角形经ETP投影后的平面坐标(x，y)之间的转换关系，利用判断点与线段的位置关系代替判断最短距离来确定地址码。在保证转换效率的同时，得到了精确的转换结果，且生成的地址码具有方向性，为球面实体的邻近查询及空间分析奠定了基础。

1 改进算法原理与步骤 1.1 QTM模型的选择

本文选取改进的QTM作为研究模型。该模型以球内切正八面体作为球面格网划分的基础，在按八面体投影划分时，能与经纬度格网中的标志性要素（如极点、赤道、主子午线等）重合，更容易确定球面上点要素的归属八面体，便于QTM格网与常用的球面经纬度格网之间的转化^[21]。

细分方法采用纬线法^[22]，即对QTM模型中的每一个球面三角形，取其三边中点作为划分点，中点连线的球面投影作为划分线，把球面三角形一分为四，依次类推，根据剖分的层次对整个球面进行近似的均匀划分，其中格网的底（顶）边用纬线代替，这样可以减小格网的变形量，提高了格网数据与地理坐标之间的转换效率，是一种更理想的QTM模型。

1.2 QTM地址编码

在QTM层次数据结构中，点的位置由QTM格元的中心点确定，随着剖分层次的增加，点位精度越来越高。任意一对经纬度坐标有且只有一个三角形地址码与之对应；反之，对于任意一个QTM地址码，就一定有一对经纬度坐标在相应的精度范围内与之对应。一个QTM位置码是由八分数（0~7）和最多30位叶结点四分数字(0~3)组成，在第k个剖分层次，三角形A的地址码表示为L_k=a₀a₁a₂a₃…a_k。其中，a₁到a_k是k个四分码，a₀是八分码，表达最初在0层次的球面正八面体剖分。由于地球被剖分为8个等球面三角形，a₀用0~7这8个数字来表达每个三角形面的地址码。对于以下递归划分的三角四叉树码a₁a₂a₃…a_k，中心三角面的ID码为0，上或下三角面的ID码为1，左三角面的ID码为2，右三角面的ID码为3。这样进行的地址码具有固定的“上下左右”方向性，有利于三角形格网的邻近查询。

1.3 改进算法原理与具体步骤

根据行列逼近法中的方法对行和列进行定义：行对应与纬线平行的一条格网离散线，列是一行内从经度为0°开始的三角形个数。算法原理如下：

1) 由地址码向经纬度转换的过程中，首先根据行列逼近算法的原理求出该地址码对应的行号i和列号j；然后应用行列号(i，j)解算出对应的经纬度(λ，φ)。

2) 由经纬度向地址码转换的过程中，首先根据待求点的纬度φ得到行号i；然后结合已知的地址码L_k-1和行号i求列号j；再将经纬度(λ，φ)转换到平面上得到平面坐标(x，y)，并利用行列号i、j以及地址码L_k-1求出0号三角形格网中两条非纬线的格边在平面上的线段方程；最后通过比较行号i的奇偶性以及平面点(x，y)与格网边线段的位置关系来确定地址码L_k。

为计算方便，可将其他八分体单元的数据首先转换至0号单元，再进行经纬度与地址码之间的转换。

1.3.1 地址码向经纬度转换

1) 根据三角形的层次K，可以求解出最大行数I：

$ I = {2^K}$

(1)

2) 对地址码L_k按一定顺序进行逐层递归求解出行号i。该过程省去了行列逼近法中借助最大行号I的递归，转换计算更加简单，效率更高。由于三角形格网存在格网方向不一致的问题，需要在层次解算中设置一个标识变量p来标识格元的方向。当p=1时，为顶三角形；当p=－1时，为底三角形。

(1) 当a_k=1时，若p=1，i=2×(i－1)+1；若p=－1，i=2×(i－1)+2。

(2) 当a_k=2或3时，若p=1，i=2×(i－1)+2；若p=－1，i=2×(i－1)+1。

(3) 当a_k=0时，若p=1，i=2×(i－1)+2；若p=－1，i=2×(i－1)+1，p=－p。

(4) 循环步骤(1)~(3)，直到求出三角形的行数i。

3) 结合三角形的行数i和剖分层次，对地址码进行逐层递归，求解格元的列数j。具体步骤沿用行列逼近法中的转换步骤，详见文献[18]。

4) 由行列号i、j求得对应3个格网点的经纬度坐标，再求其算术平均值得到格网参考点的经纬度坐标，如图 1所示。

(1) 计算格网点经度：

三角形为正向时(i为奇数)：

$ {L_T} = 90^\circ \times \frac{{\left( {j - 1} \right)/2}}{{I - i}}$

(2)

$ {L_L} = 90^\circ \times \frac{{\left( {j - 1} \right)/2}}{{I - i + 1}}$

(3)

$ {L_R} = 90^\circ \times \frac{{\left( {j - 1} \right)/2 + 1}}{{I - i + 1}}$

(4)

式中，L_T、L_L、L_R分别为三角形3个顶点的经度。

三角形为反向时(i为偶数)：

$ {L_T} = 90^\circ \times \frac{{j/2}}{{I - i + 1}}$

(5)

$ {L_L} = 90^\circ \times \frac{{j/2 - 1}}{{I - i}}$

(6)

$ {L_R} = 90^\circ \times \frac{{j/2}}{{I - i}}$

(7)

(2) 计算格网点纬度：

$ {L_U} = 90^\circ \times \frac{i}{I}$

(8)

$ {L_D} = 90^\circ \times \frac{{i - 1}}{I}$

(9)

式中，L_U、L_D分别为顶边、底边纬度。

(3) 计算格网点经纬度(λ，φ)的算术平均值：

$ \lambda = \frac{{{L_T} + {L_L} + {L_R}}}{3}$

(10)

三角形为正向时(i为奇数)：

$ \varphi = \left( {{L_U} + 2{L_D}} \right)/3$

(11)

三角形为反向时(i为偶数)：

$ \varphi = \left( {2{L_U} + {L_D}} \right)/3$

(12)

1.3.2 经纬度向地址码转换

1) 根据三角形的层次K可以求解出最大行数I。

2) 定义N为地址码中除八分码外0的个数。当N为奇数时，为正向格网；当N为偶数时，为反向格网。

$ N = \mathop {\sum\limits_{i = 1} }\limits^{k - 1} \left[ {{a_i} = 0} \right]$

(13)

式中，[a_i=0]为第i位地址码等于0的情况，若等于0，则个数加1。

3) 根据纬度坐标φ求解格网的行号i，并判断其奇偶性。

$ i = {\rm{int}}\left[ {\frac{{\varphi \times I}}{{90}}} \right] + 1$

(14)

当格网为正向(N为奇数)且i为偶数时，或当格网为反向(N为偶数)且i为奇数时，a_k=1；否则进行下一步计算。

4) 将经纬度坐标(λ，φ)转换为平面坐标(x，y)：

$ x = \frac{I}{{\rm{ \mathit{ π} }}}\left[ {\varphi + 2\lambda \left( {1 - \frac{2}{{\rm{ \mathit{ π} }}}\varphi } \right)} \right]$

(15)

$ y = \frac{{\sqrt 3 I}}{{\rm{ \mathit{ π} }}}\varphi $

(16)

5) 将前一层地址码L_k-1加0，得到第k层对应四元格网中0号格网单元的地址码code0，由code0与行号i（参考前文转换过程）得到对应列号j。

6) 由行列号i、j经数学运算得到0号格网底边（顶边）两格网点的平面坐标分别为P₁(x，y)、P₂(x，y)。

7) 在平面坐标系中，0号格网中两格网边的直线方程为L_l、L_r（以正向格网为例）：

$ {L_l}:y = - \sqrt 3 x + {b_1}$

(17)

$ {L_r}:y = \sqrt 3 x + {b_2}$

(18)

因为P₁、P₂分别为L_l、L_r上的点，可利用式(17)和式(18)求得b₁、b₂。

8) 根据待求点的x坐标判断所属格网，如图 2所示。

用P₁.x、P₂.x分别表示P₁、P₂点在平面坐标系中的x坐标，以下仅以正向三角格网为例，判断所属格网的过程，倒向三角格网的判断原理相同。

(1) 当x < P₁.x时，a_k=2。

(2) 当x > P₂.x时，a_k=3。

(3) 当P₁.x < x < (P₁.x+P₂.x)/2时，将x代入式(17)得到y₁，若y > y₁，则a_k=0；若y < y₁，则a_k=2。

(4) 当(P₁.x+P₂.x)/2 < x < P₂.x时，将x代入式(18)得到y₂，若y > y₂，则a_k=0；若y < y₂，则a_k=3。

9) 重复步骤1)~8)，循环次数由所需地址码位数决定。

改进后的算法利用判断点与线段的位置来代替ETP投影法中判断最短距离的过程，在计算过程中只有加、减、乘、除等简单的矢量算术运算，在保证准确转换的基础上大大提高了转换效率，且生成的地址码依然具有方向性，能有效进行邻近查询。

2 实验对比分析

为了比较不同算法间的转换效果，本文选取不同层次（12、14、16、19、21层，分别对应1:400万、1:100万、1:25万、1:5万、1:1万比例尺）、不同数量的坐标点（1万、10万、100万、1 000万），依次利用5种算法（ETP投影法、ZOT投影法、行列逼近法、三向转换法、改进算法）进行地址码与经纬度之间的相互转换。采用的实验设备为ASUS X550V笔记本；硬件配置为Inter(R) Core(TM) i5-3230M CPU @ 2.60 GHz，4 GB RAM；操作系统为Windows 10专业版；开发工具为Visual Studio 2010。

2.1 计算效率对比分析

5种算法效率的实验结果如表 1所示。通过表 1可以看出，改进算法的转换效率达到该类算法的最优级别，与ZOT投影法、行列逼近法、三向转换法处于同一量级。改进算法的时间消耗为ETP投影法的10.1%（地址码→经纬度，见图 3(a)）和10.4%（经纬度→地址码，见图 3(b)）。

2.2 算法精度分析

在转换精度方面，保证了转换效率以及方向性的改进算法同样可以得到准确的转换结果。由经纬度到地址码的转换过程中，行列逼近法较改进算法、ETP投影法和ZOT投影法多出半个格网误差，由地址码转换所得的经纬度坐标没有固定规则，不具参考性。三向转换算法在确定三轴坐标的过程中，将所有格网边均作为小圆弧处理，由于球面为非欧空间，当格网边均以小圆弧处理时，不可避免地存在边界“交叉地带”，从而造成转换误差，且影响范围随着剖分层次的增加而增大。改进算法中，格网的参考点是通过格网点坐标计算所得，可根据不同应用需求改变参考点的计算公式，灵活选择参考点。

3 结语

本文针对已有算法存在的问题，提出了一种改进的QTM地址码与经纬度转换算法，并给出了详细的算法原理及转换步骤，通过实验对比得出如下结论：

1）改进算法在行列逼近的基础上加入判断点与线段位置关系的操作，在不丢失效率的基础上得到精确的转换结果；解决了行列逼近法和三向转换法中无法保证转换精度的问题；且改进算法可根据应用需求，灵活确定转换所得的参考点位置，更具普遍适用性。

2）与ETP投影法相比，改进算法的计算过程仅涉及简单的算术运算，转换效率大大提高，转换时间消耗与ZOT投影法和行列逼近法处于同一量级，其中经纬度→地址码转换过程仅为ETP投影法时间消耗的10.4%，地址码→经纬度转换过程仅为ETP投影法时间消耗的10.1%。

3）改进算法得到的地址码同样具有方向性，利于计算机的邻近搜索等操作。