本發(fā)明涉及地理信息系統(tǒng)技術，具體涉及一種球面三角形離散格網(wǎng)編碼向地理經(jīng)緯度坐標的快速轉(zhuǎn)換方法。

背景技術：

球面三角形離散格網(wǎng)是一種經(jīng)典的全球離散格網(wǎng)系統(tǒng)。它以三角形為基本格網(wǎng)單元，基于正八面體或正二十面體對球面空間進行同構(gòu)離散化處理，形成具有層次遞歸性、單元形態(tài)均勻性，無縫無疊的多尺度全球格網(wǎng)系統(tǒng)(tong，2013)。經(jīng)過眾多學者長期的研究，球面三角形離散格網(wǎng)不僅成為了全球多分辨率空間數(shù)據(jù)建模的有效工具，還被廣泛應用于空間數(shù)據(jù)索引與可視化(goodchild，1992；feket，1990)、地圖綜合模型(dutton,1996)、全球?qū)Ш侥Ｐ?lee，2000)、全球環(huán)境監(jiān)測(dutton，1999)、全球氣候模擬(sadourny，1968；heikes，1995；1996；thuburn，1997)等領域。

編碼作為球面三角形離散格網(wǎng)單元的唯一數(shù)字標識，既表達了空間位置信息，又蘊含了空間尺度信息，為處理全球多分辨率海量空間數(shù)據(jù)提供了可能，是格網(wǎng)系統(tǒng)的重要組成部分(dutton，1989，2000；fekete，1990；goodchild，1992；otoo，1993；lee，2000；bartholdi，2001)。然而，現(xiàn)有的大部分空間數(shù)據(jù)是以地理坐標為主，許多傳統(tǒng)空間信息系統(tǒng)也是以地理坐標為基礎。因此，如何準確、高效的實現(xiàn)球面三角形離散格網(wǎng)編碼與地理坐標的雙向轉(zhuǎn)換是球面三角形離散格網(wǎng)研究中的核心問題之一(dutton，1999)。

由于球面是一個各向異性的二維緊致流形空間，在拓撲性質(zhì)上與平面空間不同胚，屬于不可展曲面(kimerling，1999)。因此，球面空間無法與平面空間一樣用相同的格網(wǎng)連續(xù)鋪砌，難以借助格網(wǎng)原點與目標點的距離和角度關系直接推算出目標點對應的格網(wǎng)編碼，通常采用遞歸逼近的方式逐層建立地理坐標和格網(wǎng)編碼的對應關系。在此方面，相關學者基于正八面體球面三角形格網(wǎng)，先后提出了zot投影(dutton，1991)、etp投影(goodchild，1992)、sqc(otoo，1993)、rca(zhao，2006)、三向互化(tong，2006)、eatp投影(sun，2007)等方法。這些算法充分利用了基于正八面體的球面三角形格網(wǎng)與地理坐標具有簡易轉(zhuǎn)換的特性(sahr，2003)，從不同角度簡化了遞歸逼近時的算法實現(xiàn)難度，提高了雙向轉(zhuǎn)換的效率。

然而，這些算法的本質(zhì)仍是基于遞歸逼近的方法逐層建立地理坐標與格網(wǎng)編碼的映射關系，算法時間復雜度與格網(wǎng)剖分層次呈線性關系。隨著格網(wǎng)剖分層次的提高，算法效率明顯下降，不利于高分辨率空間數(shù)據(jù)的集成。同時，這些改進方法僅適用于格網(wǎng)結(jié)構(gòu)簡單但幾何性質(zhì)較差的正八面體球面三角形格網(wǎng)，難以應用在格網(wǎng)結(jié)構(gòu)復雜但幾何性質(zhì)良好的正二十面體球面三角形格網(wǎng)上，不具有通用性。

通過分析可以發(fā)現(xiàn)，盡管現(xiàn)有算法充分考慮到球面空間的特殊性對球面三角形離散格網(wǎng)系統(tǒng)中格網(wǎng)單元編碼與地理坐標雙向轉(zhuǎn)換的影響，但卻忽略了這一特殊性質(zhì)在格網(wǎng)剖分層次提高過程中的變化規(guī)律。球面三角形離散格網(wǎng)的幾何性質(zhì)并不是固定不變的。球面三角形離散格網(wǎng)在最初的剖分層次中格網(wǎng)單元幾何性質(zhì)確實差異較大，但隨著剖分層次的提高，則幾何性質(zhì)趨于相同(white，1998)，局部區(qū)域的幾何性質(zhì)已近似于平面。若能根據(jù)格網(wǎng)幾何性質(zhì)的變化規(guī)律，尋找合適的控制剖分層次作為閾值，在低于該閾值的剖分層次中采用可以確保精度但效率較低的遞歸逼近轉(zhuǎn)換方法，在高于該閾值的剖分層次采用類似平面格網(wǎng)的直接映射方法，則有可能構(gòu)建一種兼顧精度與效率的新型混合式轉(zhuǎn)換算法。

技術實現(xiàn)要素：

發(fā)明目的：本發(fā)明的目的在于解決現(xiàn)有技術中存在的不足，提供一種球面三角形離散格網(wǎng)編碼向地理經(jīng)緯度坐標的快速轉(zhuǎn)換方法。

技術方案：本發(fā)明一種球面三角形離散格網(wǎng)編碼向地理經(jīng)緯度坐標的快速轉(zhuǎn)換方法，依次包括以下步驟：

(1)輸入格網(wǎng)編碼；

(2)根據(jù)編碼的剖分層次l確定控制剖分層次lc；

(3)在低于控制剖分層次的階段，采用遞歸逼近方法確定格網(wǎng)編碼在控制剖分層次上所對應的控制三角形單元v0v1v2及目標編碼格元在控制三角形v0v1v2內(nèi)的行列號(x，y)；

(4)當轉(zhuǎn)換到控制剖分層次lc時，采用直接映射法將控制三角形v0v1v2采用等分方式直接劃分成格網(wǎng)，而后根據(jù)控制三角形內(nèi)的行列號(x，y)找到目標三角形格元；

(5)計算目標三角形格元中心點坐標，將其從直角坐標轉(zhuǎn)換到經(jīng)緯度坐標，輸出完成轉(zhuǎn)換。

進一步的，所述步驟(3)中在低于控制剖分層次的階段時，采用遞歸逼近法確定格網(wǎng)編碼在控制剖分層次上所對應的控制三角形單元v0v1v2目標編碼格元在控制三角形v0v1v2內(nèi)的行列號(x，y)，具體過程為：

(3.1)根據(jù)格網(wǎng)編碼獲取對應初始三角形單元(例如正八面體球面三角形離散格網(wǎng)的初始八個球面三角形，或正二十面體球面三角形離散格網(wǎng)的球面三角形)，用其頂點初始化控制三角形v0v1v2；

(3.2)初始化格網(wǎng)剖分計數(shù)l＝0，根據(jù)各種編碼與行列號的對應關系(例如二進制sqc編碼，抽出其morton碼的奇數(shù)位來得到的新二進制數(shù)就是列號y，morton碼的偶數(shù)位組成的二進制數(shù)加上方位碼就是行號x)，將格網(wǎng)編碼轉(zhuǎn)換為行列號(x，y)；

(3.3)若當前控制三角形的層次l等于控制剖分層次lc，輸出控制三角形單元v0v1v2，及剩余(x，y)，遞歸算法結(jié)束；否則，就進入步驟(3.4)；此處剩余(x，y)是指控制三角形內(nèi)的行號和列號；

(3.4)計算控制三角形單元三條弧邊的中點m0、m1、m2，通過x、y判斷下一層次控制三角形在當前控制三角形中的位置，若y<2^l-l-1，則為當前控制三角形的上部子三角形，進入步驟(3.5)；若x<2y-2^l-l，則為當前控制三角形的左部子三角形，令y＝y(tǒng)–2^l-l-1，進入步驟(3.5)；若y<2^l-l-1，則為當前控制三角形的中部子三角形，令x＝x–2×y+2^l-l-1，y＝2^l-l-1-y-1，進入步驟(3.5)；否則，為當前控制三角形的上部子三角形，令x＝x-2^l-l，y＝y(tǒng)-2^l-l-1，進入步驟(3.5)；

(3.5)更新控制三角形單元的頂點坐標v0、v1、v2，剖分計數(shù)l＝l+1；返回步驟(3.3)。

進一步的，所述步驟(4)中步驟(2)通過遞歸逼近方法獲得控制剖分層次上的控制三角形單元v0v1v2后，根據(jù)剩余的剖分次數(shù)n＝l-lc,采用直接映射方法得到目標剖分層次上的三角形單元，參照平面柵格，將控制三角形v0v1v2采用等分方式直接劃分成格網(wǎng)，而后根據(jù)控制三角形內(nèi)的行列號(x，y)找到目標三角形格元，具體流程如下：

(4.1)將v0v1、v0v2等分成2ⁿ段，在v0v1上分別取第y和y+1個等分點e1、e2，在v0v2上分別取第y，y+1個等分點e3、e4；

(4.2)將e1e3等分成y段，e2e4等分成y+1段，而后在e1e3上分別取第個等分點d1、d2，在e2e4上分別取第個等分點d3、d4；

(4.3)若x為偶數(shù)，則目標層次上的三角形三個頂點為d1、d3、d4；若為奇數(shù)，則三個頂點為d4、d1、d2；

進一步的，基于步驟(4.3)所得的三個頂點計算三角形中心點坐標，并將其從直角坐標轉(zhuǎn)換到經(jīng)緯度坐標，完成轉(zhuǎn)換。

有益效果：依據(jù)三角形格網(wǎng)的幾何特征，通過遞歸逼近方法獲得控制剖分層次上的控制三角形單元v0v1v2后，根據(jù)剩余的剖分次數(shù)n＝l-lc，采用直接映射方法得到目標剖分層次上的三角形單元，參照平面柵格，將控制三角形v0v1v采用等分方式直接劃分成格網(wǎng)，而后根據(jù)控制三角形內(nèi)的行列號(x，y)找到目標三角形格元。

附圖說明

圖1為本發(fā)明的整體流程示意圖；

圖2為本發(fā)明中基于遞歸逼近方法的格網(wǎng)編碼向地理坐標的轉(zhuǎn)換流程示意圖；

圖3為實施例中三角形格網(wǎng)的sqc編碼與行列號對應關系示意圖；

圖4為實施例中子三角形位置分布示意圖；

圖5為實施例中控制三角形單元內(nèi)的等分格網(wǎng)。

具體實施方式

下面對本發(fā)明技術方案進行詳細說明，但是本發(fā)明的保護范圍不局限于所述實施例。

如圖1所示，本發(fā)明的一種球面三角形離散格網(wǎng)編碼向地理經(jīng)緯度坐標的快速轉(zhuǎn)換方法，依次包括以下步驟：

(1)輸入格網(wǎng)編碼；

(2)根據(jù)編碼的剖分層次l確定控制剖分層次lc；

(3)在低于控制剖分層次的階段，采用遞歸逼近方法確定格網(wǎng)編碼在控制剖分層次上所對應的控制三角形單元v0v1v2及目標編碼格元在控制三角形v0v1v2內(nèi)的行列號(x，y)；

(4)當轉(zhuǎn)換到控制剖分層次lc時，采用直接映射法將控制三角形v0v1v2采用等分方式直接劃分成格網(wǎng)，而后根據(jù)控制三角形內(nèi)的行列號(x，y)找到目標三角形格元；

(5)計算目標三角形格元中心點坐標，將其從直角坐標轉(zhuǎn)換到經(jīng)緯度坐標，輸出完成轉(zhuǎn)換。

所述步驟(3)中在低于控制剖分層次的階段時，采用遞歸逼近法確定格網(wǎng)編碼在控制剖分層次上所對應的控制三角形單元v0v1v2目標編碼格元在控制三角形v0v1v2內(nèi)的行列號(x，y)，具體過程為：

(3.1)根據(jù)格網(wǎng)編碼獲取對應初始三角形單元(正八面體球面三角形離散格網(wǎng)的初始八個球面三角形，或正二十面體球面三角形離散格網(wǎng)的球面三角形)，用其頂點初始化控制三角形v0v1v2；

(3.2)初始化格網(wǎng)剖分計數(shù)l＝0，如圖2所示，sqc編碼與行列號的對應關系，將格網(wǎng)編碼轉(zhuǎn)換為行列號(x，y)；

(3.3)若當前控制三角形的層次l等于控制剖分層次lc，輸出控制三角形單元v0v1v2，及剩余(x，y)，遞歸算法結(jié)束；否則，就進入步驟(3.4)；此處剩余(x，y)是指控制三角形內(nèi)的行號和列號；

(3.4)計算控制三角形單元三條弧邊的中點m0、m1、m2，通過x、y判斷下一層次控制三角形在當前控制三角形中的位置，若y<2^l-l-1，則為當前控制三角形的上部子三角形，進入步驟(3.5)；若x<2y-2^l-l，則為當前控制三角形的左部子三角形，令y＝y(tǒng)–2^l-l-1，進入步驟(3.5)；若y<2^l-l-1，則為當前控制三角形的中部子三角形，令x＝x–2×y+2^l-l-1，y＝2^l-l-1-y-1，進入步驟(3.5)；否則，為當前控制三角形的上部子三角形，令x＝x-2^l-l，y＝y(tǒng)-2^l-l-1，進入步驟(3.5)；

(3.5)更新控制三角形單元的頂點坐標v0、v1、v2，剖分計數(shù)l＝l+1；返回步驟(3.3)。

所述步驟(4)中步驟(2)通過遞歸逼近方法獲得控制剖分層次上的控制三角形單元v0v1v2后，根據(jù)剩余的剖分次數(shù)n＝l-lc,采用直接映射方法得到目標剖分層次上的三角形單元，參照平面柵格，將控制三角形v0v1v2采用等分方式直接劃分成格網(wǎng)，如圖4所示，而后根據(jù)控制三角形內(nèi)的行列號(x，y)找到目標三角形格元，具體流程如下：

(4.1)將v0v1、v0v2等分成2ⁿ段，在v0v1上分別取第y和y+1個等分點e1、e2，在v0v2上分別取第y，y+1個等分點e3、e4；

(4.2)將e1e3等分成y段，e2e4等分成y+1段，而后在e1e3上分別取第個等分點d1、d2，在e2e4上分別取第個等分點d3、d4；

(4.3)若x為偶數(shù)，則目標層次上的三角形三個頂點為d1、d3、d4；若為奇數(shù)，則三個頂點為d4、d1、d2；

最后，基于步驟(4.3)所得的三個頂點計算三角形中心點坐標，并將其從直角坐標轉(zhuǎn)換到經(jīng)緯度坐標，完成轉(zhuǎn)換。

實施例1：

本實施例中，以南京師范大學仙林校區(qū)地理科學學院所在的正二十面體三角形離散格網(wǎng)22層格元的sqc編碼為例，使用本發(fā)明算法將其轉(zhuǎn)換成地理坐標，該sqc編碼為(110111010101101000010100100101011001111000000)。

本實施例的轉(zhuǎn)換過程分為兩個重要部分，控制剖分層次前的遞歸逼近和控制剖分層次后的直接映射。

一、遞歸逼近過程步驟如下：

步驟1：根據(jù)格網(wǎng)編碼獲取對應編號為2的初始三角形單元，用其頂點初始化控制三角形v0v1v2；

步驟2：初始化格網(wǎng)剖分計數(shù)l＝0，如圖2所示，sqc編碼與行列號的對應關系，將格網(wǎng)sqc編碼轉(zhuǎn)換為行列號(5342296，4135384)；

步驟3：若當前控制三角形的層次l等于控制剖分層次lc，輸出控制三角形單元v0v1v2，及剩余行列號(x，y)(即控制三角形內(nèi)的行列號)，遞歸算法結(jié)束；否則，就進入下一步；

步驟4：計算控制三角形單元三條弧邊的中點m0、m1、m2，通過x、y判斷下一層次控制三角形在當前控制三角形中的位置，若y<2^l-l-1，則為當前控制三角形的上部子三角形，進入步驟5；若x<2y-2^l-l，則為當前控制三角形的左部子三角形，令y＝y(tǒng)–2^l-l-1，進入步驟5；若y<2^l-l-1，則為當前控制三角形的中部子三角形，令x＝x–2×y+2^l-l-1，y＝2^l-l-1-y-1，進入步驟5；否則，為當前控制三角形的上部子三角形，令x＝x-2^l-l，y＝y(tǒng)-2^l-l-1，進入步驟5。

步驟5：更新控制三角形單元的頂點坐標v0、v1、v2，剖分計數(shù)l＝l+1；返回步驟3。

步驟3到步驟5是一個遞歸逼近過程，每次逼近行列號(x,y)的值如下：(5342296，4135384)、(1147992，2038232)、(1147992，989656)、(99416，465368)、(99416，203224)、(99416，72152)、(86183，58919)、(20647，26151)。

同時得到最終的控制三角形單元v0v1v2，其三個頂點的空間直角坐標分別為

v0(-2627.6295552955821，4766.8453899067490，3366.4147024485728)、

v1(-2599.5515195261419，4748.9227890212578，3413.2190438417356)、

v2(-2649.1708091961154，4724.8164347521324，3408.5192799219058)。

二、直接映射過程步驟如下：

步驟1：將v0v1、v0v2等分成2¹⁵段，在v0v1上分別取第26151，26152個等分點e1、e2，在v0v2上分別取第26151，26152個等分點e3、e4。這幾點的坐標為：

e1(-2605.2383667798254,4752.5728418204199,3403.7897166907010)、

e2(-2605.2375079911003,4752.5722913790933,3403.7911425603065)、

e3(-2644.8417093627831,4733.3408211258838,3400.0436767015644)、

e4(-2644.8423644033683,4733.3395342855165,3400.0449586184859)。

步驟2：將e1e3等分成26151段，e2e4等分成26152段，而后在e1e3上分別取第10323，10324個等分點d1、d2，在e2e4上分別取第10323，10324個等分點d3、d4。這幾點的坐標為：

d1(-2620.8865472413713,4745.0080987270630,3402.3303557273480)、

d2(-2620.8880621629605,4745.0073642110137,3402.3302131326013)、

d3(-2620.8872040829283,4745.0068150879670,3402.3316399556224)、

d4(-2620.8856891612695,4745.0075496038899,3402.3317825502772)。

步驟3：因為x為奇數(shù)，則目標層次上的三角形三個頂點為d4、d1、d2；

最后，基于上述步驟所得的三個頂點計算三角形中心點坐標d(-2620.8872711624540,4745.0074260087431,3402.3307362719015)，并將其從直角坐標轉(zhuǎn)換到經(jīng)緯度坐標，完成轉(zhuǎn)換，最終所得經(jīng)緯度坐標為：p(118.91390283576838°e,32.114583311947690°n)。

通過上述實施例可以看出，本發(fā)明利用球面三角形離散格網(wǎng)的幾何性質(zhì)的變化規(guī)律，這一規(guī)律表明球面三角形離散格網(wǎng)的幾何性質(zhì)是不斷變化的，隨著剖分層次的提高不斷趨于近似。在低剖分層次中，由于格網(wǎng)單元幾何性質(zhì)差異較大，且格網(wǎng)曲率與平面差異明顯，因此采用層次遞歸逼近轉(zhuǎn)換法以確保編碼與解碼的正確性；在高剖分層次中，由于格網(wǎng)單元幾何性質(zhì)近似相同，格網(wǎng)曲率近乎于平面，因此采用類似平面格網(wǎng)的直接映射方法，在確保轉(zhuǎn)換精度的前提下，提高轉(zhuǎn)換的效率。