耿曉暉,程承旗,宋樹華,李大鵬

(北京大學(xué)遙感與地理信息系統(tǒng)研究所,北京 100871）

基于地圖方里網(wǎng)的全球剖分系統(tǒng)

耿曉暉,程承旗,宋樹華,李大鵬

(北京大學(xué)遙感與地理信息系統(tǒng)研究所,北京 100871）

針對現(xiàn)有剖分模型的不足,提出了一種基于地圖方里網(wǎng)的全球剖分系統(tǒng),有效避免了傳統(tǒng)經(jīng)緯度格網(wǎng)模型在高緯度地區(qū)的形狀退化和正多面體格網(wǎng)模型的面片形狀不規(guī)則問題。制定了相應(yīng)的編碼,實現(xiàn)了面片編碼與傳統(tǒng)地理坐標(biāo)之間的轉(zhuǎn)換和鄰接關(guān)系的計算,最后對方里網(wǎng)在地球橢球面上的變形規(guī)律進行了研究。研究結(jié)果表明:基于地圖方里網(wǎng)的剖分系統(tǒng)不僅具有科學(xué)的數(shù)據(jù)組織形式,而且面片單元的變形面積小、變形規(guī)律穩(wěn)定。

全球剖分系統(tǒng);方里網(wǎng);高斯克呂格投影;Mo rton碼

0 引言

隨著空間技術(shù)和對地觀測技術(shù)的飛速發(fā)展,人們能夠更加方便地獲取有關(guān)地球及其各種資源、環(huán)境和社會現(xiàn)象的多分辨率的、海量的、實時的對地觀測數(shù)據(jù)[1]。面對海量空間數(shù)據(jù)管理問題,傳統(tǒng)的空間信息表達、組織、管理和發(fā)布方式已不能滿足全球數(shù)據(jù)管理的要求。因此,需要構(gòu)建一個新的基于全球的、多尺度的、融合了空間索引機制的、無縫的、開放的層次性空間數(shù)據(jù)管理框架[2-4]。基于各種格網(wǎng)模型的全球剖分系統(tǒng)成為解決這一問題的新思路,現(xiàn)有的剖分系統(tǒng)主要有經(jīng)緯度格網(wǎng)剖分系統(tǒng)、正多面體格網(wǎng)剖分系統(tǒng)和自適應(yīng)格網(wǎng)剖分系統(tǒng)[1]。鄧雪清等采用經(jīng)緯度格網(wǎng)模型剖分的方法[5,6],該方法計算簡單,保持與原有地理數(shù)據(jù)格式的兼容性,但在高緯度地區(qū)存在明顯的形狀退化。趙學(xué)勝等采用正多面體格網(wǎng)模型剖分的方法[7,8],該方法避免了高緯度地區(qū)的形狀退化,但面片多為三角形、菱形和六邊形,不符合計算機的表達習(xí)慣,往往運算復(fù)雜。自適應(yīng)格網(wǎng)模型的面片形狀取決于采樣點的分布,缺乏幾何規(guī)則性。

將方里網(wǎng)作為剖分模型,能夠有效地避免上述模型存在的缺點:方里網(wǎng)不但在平面內(nèi)形狀統(tǒng)一、量測方便,而且根據(jù)高斯投影的變形規(guī)律,緯度越大并行越小,所以不存在兩極退化;同時高斯投影屬于等角投影,即相互垂直的兩條直線投影之后夾角仍為90°,保持了面片角點的正方形特性。本文以方里網(wǎng)為基礎(chǔ),構(gòu)建一種基于地圖方里網(wǎng)的全球剖分系統(tǒng)。

1 方里網(wǎng)作為剖分模型的難點和解決方法

方里網(wǎng)是一種建立在某種地圖投影基礎(chǔ)上的格網(wǎng)系統(tǒng),將制圖區(qū)域按平面坐標(biāo)或按經(jīng)緯度劃分為格網(wǎng),以格網(wǎng)為單位描述或表達其中的屬性分類、統(tǒng)計分級以及變化參數(shù),即在二維上表達動態(tài)時空變化的規(guī)律[9]。方里網(wǎng)作為全球剖分模型存在的難點主要是跨帶問題。如圖1所示,在平面坐標(biāo)系上投影帶之間存在縫隙,縫隙間的格網(wǎng)屬于外部面片,如果利用高斯公式直接反算至地球橢球面,會造成鄰帶方里網(wǎng)之間的重疊、壓蓋,如果通過邊界經(jīng)線切割后直接去除,會破壞投影帶在平面上的連續(xù)性,增加統(tǒng)一編碼的難度。

圖1 方里網(wǎng)面片的跨帶Fig.1 Kilo-grids step across the projection zone

本文采用的方法是:首先在平面內(nèi)對所有的面片進行統(tǒng)一編碼,然后判斷面片與投影帶的關(guān)系,并以屬性的方式進行標(biāo)識。方里網(wǎng)面片分為3種:1）內(nèi)部面片。面片的4個角點在投影帶內(nèi)部,可以通過高斯反解公式直接計算至球面[9]。2）外部面片。面片4個角點在投影帶外部,不做至球面的計算,只存儲一個編碼,保證平面計算的連續(xù)性,所以也可以認(rèn)為是虛擬面片。3）邊界面片。面片的部分角點在投影帶內(nèi),需要投影帶的切割,只將在投影帶內(nèi)的部分反算至球面。最后將面片的邊界屬性存儲至面片編碼。在計算過程中首先進行邊界屬性判斷,根據(jù)面片的邊界屬性分別做不同處理,保證在平面內(nèi)編碼的統(tǒng)一,同時避免了將外部的方里網(wǎng)反算至地球球面造成方里網(wǎng)重疊。

2 基于方里網(wǎng)的球面剖分系統(tǒng)

2.1 剖分系統(tǒng)的建立和編碼

基于上述解決方案,筆者提出構(gòu)建基于地圖方里網(wǎng)的全球剖分系統(tǒng)的總體思路:以1 km×1 km的方里網(wǎng)為基礎(chǔ),以四叉樹的形式向上綜合、向下等分,編碼的方式采用Mo rton碼。

Morton碼是一種線性的四叉樹編碼方式,其基本思想是:不需記錄中間結(jié)點和使用指針,僅記錄葉結(jié)點,并用地址碼表示葉結(jié)點的位置,編碼規(guī)則如圖2,編碼與行列數(shù)間的轉(zhuǎn)換表示為[10]:

式中:i,j表示面片行列,m(i,j）表示面片的Morton碼:int(）為向下取整運算,mod為求余函數(shù)。

圖2 Morton碼的編碼規(guī)則Fig.2 Coding rules of Morton

剖分流程如下:1）投影帶擴展:經(jīng)過6°分帶后,投影帶中央經(jīng)線的長度為2×10 002 km,在赤道的寬度為2×334 km,為了能夠以四叉樹的形式組織面片,投影帶的長寬必須是2nkm,所以首先將投影帶的長度擴展為2×10 240 km,寬度為2×512 km。2）初級面片劃分:首先通過投影坐標(biāo)軸將擴展后的投影帶劃分為4個象限,然后將每個象限等分為20行,這樣整個區(qū)域被劃分為4×20個512 km×512 km的面片,其編碼方式為:象限號(quadrant）、行號(row）。3）次級面片劃分:以初級面片為基礎(chǔ)以四叉樹的形式向下等分,得到4個大小為256 km× 256 km的面片,作為第二級面片,編碼方式采用Morton碼。依此類推可以得到大小為128 km× 128 km、64 km×64 km、…的第二級、第三級、…、第n級面片(圖3）。4）邊界屬性判斷:通過高斯反解公式[9]計算面片的4個角點坐標(biāo)經(jīng)度值,如果全部在兩條邊界經(jīng)線的經(jīng)度范圍之間,為內(nèi)部面片,全部超出該范圍為外部面片,其余為邊界面片。

圖3 基于方里網(wǎng)剖分系統(tǒng)的面片劃分規(guī)則和編碼Fig.3 Division and coding rules of the subdivision system based on Kilo-grids

方里網(wǎng)剖分面片編碼包括5項(圖4）:1）面片級別(level）:定義最高的面片為20級(第20級為米級）,用5位二進制數(shù)表示。2）投影帶號(zone）:全球共60帶,用6位二進制數(shù)表示。3）初級面片編碼(prim ary）:包括象限號(quadrant）、行號(row）。其中,象限數(shù)為4,可用2位二進制數(shù)表示,行數(shù)最大值為20,則可用5位二進制數(shù)表示,所以初級面片總共可以用7位二進制數(shù)表示。4）基于四叉樹的Morton碼:由1個變長碼組成,面片級別每增大1級, Morton增加2位,總共可以用2×level位表示。5）邊界屬性(border）:共有內(nèi)部、外部和邊界3種面片屬性,可以由2位二進制數(shù)表示。

圖4 基于方里網(wǎng)剖分面片編碼Fig.4 The code of a Kilo-grid tile

2.2 面片編碼與經(jīng)緯度之間的轉(zhuǎn)換

方里網(wǎng)面片的編碼不僅可以標(biāo)識面片,而且包含足夠的空間定位信息,可以方便地與經(jīng)緯度之間進行換算,具體換算關(guān)系如下。

通過編碼計算經(jīng)面片中心點的經(jīng)緯度坐標(biāo):

(1）初級面片編碼(primary）分為象限號(quadrant）、行號(row）,通過象限號確定高斯坐標(biāo)的正負(fù)號,計算初級面片的左下角點高斯坐標(biāo)(X0, Y0）,X0=512×row,Y0=0;

(2）用公式(2）將Morton碼轉(zhuǎn)換為面片的行列號(i,j）;

(3）計算面片中心點的高斯坐標(biāo)(X,Y）;

(4）運用高斯坐標(biāo)反算公式計算該點的經(jīng)緯度坐標(biāo)(B,L）;

(5）算法結(jié)束。

通過經(jīng)緯度求面片編碼:

(1）確定投影帶號:zone=L/6;

(2）通過高斯正解公式計算該點的高斯坐標(biāo)(X, Y）;確定初級面片:通過X、Y的正負(fù)符號確定象限數(shù),計算行數(shù),row=X/512,得到初級面片編碼;

(3）將初級面片四等分劃分為4個面片,求出相應(yīng)的編碼和4個角點的坐標(biāo),記下包含(X,Y）的面片編碼,記錄其Mo rton碼;

(4）對第 3步做level次循環(huán)即得到完整的Morton碼;

(5）判斷面片與邊界的關(guān)系,確定面片的邊界屬性(border）;

(6）算法結(jié)束。

2.3 基于面片編碼的鄰接算法

以投影帶中的第一象限為例,算法如下:

(1）通過Morton碼計算面片在投影帶中的行列數(shù)(i,j）;

(2）計算相鄰面片的行列數(shù):左面片(i-1,j）,右面片(i+1,j）,上面片(i,j+1）,下面片(i,j-1）;

(3）通過公式(1）計算鄰接面片編碼:左面片m (i-1,j）,右面片m(i+1,j）,上面片m(i,j+1）,下面片m(i,j-1）;

(4）判斷面片邊界屬性,組合各編碼分項,得到完整編碼;

(5）算法結(jié)束。

3 方里網(wǎng)面片的變形研究

方里網(wǎng)在地球橢球面上的面積變形可以通過高斯克呂格投影的長度變形公式進行計算,公式推導(dǎo)過程為:高斯克呂格投影后的長度變形公式取前兩項,即:

其中:B為緯度,l為相對經(jīng)差,η=e′cosB,e′為第二偏心率[9]。

將長度變形公式平方之后得到面積變形公式,取前兩項為:

以上是橢球面至平面的面積變形,取倒數(shù)得到平面至橢球面的面積變形:

橢球面上的微分梯形計算公式為:

式中:e為第一偏心率,a為地球橢球的長半軸[9]。

這樣方里網(wǎng)在地球橢球面上的面積變形可以通過積分計算:

其中:(B0,l0）、(B1,l1）為方里網(wǎng)投影到地球球面上的地理坐標(biāo)范圍。

定義方里網(wǎng)面積變形系數(shù)Ф:

以1 km×1 km、2 km×2 km面片為例,參數(shù)采用2000國家大地坐標(biāo)系的地球橢球參數(shù),在相應(yīng)位置的方里網(wǎng)面片的變形系數(shù)如表1、表2。

表1 1 km×1 km面片投影后的面積變形Table 1 Area distortion of the 1 km×1 km Kilo-grid tiles due to map projection

表2 2 km×2 km面片投影后的面積變形Table 2 Area distortion of the 2 km×2 km Kilo-grid tiles due tomap projection

從表1、表2得出如下結(jié)論:1）方里網(wǎng)面積變形系數(shù)全部小于等于0,即方里網(wǎng)投影到橢球面上之后面積變小;2）變形系數(shù)較小,以2 km×2 km的面片為例,最大變形數(shù)量級為10-3;3）對于同一大小的方里網(wǎng)面片,面片變形隨著相對中央經(jīng)線的經(jīng)差增加而增大,隨著緯度的增加而減小;4）方里網(wǎng)最大的變形面片在緯度為0°、經(jīng)差為3°處,所以該系統(tǒng)不適合低緯度地區(qū)的地理信息表達。

4 結(jié)語

本文通過方里網(wǎng)剖分模型,構(gòu)建了一種全新的剖分系統(tǒng),以四叉樹為基礎(chǔ),制定了相應(yīng)的編碼,研究了面片的變形。結(jié)論如下:1）編碼高效。采用變長Morton碼,只對葉結(jié)點編碼,隱含了葉結(jié)點的位置和大小,節(jié)省了存儲空間;編碼中含有面片的內(nèi)外屬性,為下一步的空間計算奠定了基礎(chǔ)。2）具有較好的量測性。方里網(wǎng)具有天然的距離和面積量測屬性,符合人們在平面的空間認(rèn)知習(xí)慣,適用于大、中比例尺空間數(shù)據(jù)的表達。3）具有等角性質(zhì)。高斯投影屬于等角投影,即角度變形為0,方里網(wǎng)投影之后,4個角仍然為直角,最大限度地保證了面片的正方形結(jié)構(gòu)。4）變形穩(wěn)定。存在面積和長度變形,符合高斯投影的變形規(guī)律。

方里網(wǎng)面片在高緯度地區(qū)變形較小,與傳統(tǒng)的經(jīng)緯度模型優(yōu)勢互補,因此構(gòu)建兩者結(jié)合的混合模型,并研究兩個系統(tǒng)之間編碼的相互轉(zhuǎn)換方法,將是下一步研究的重點。

[1] 冉令輝.全球空間信息剖分編碼模型研究[D].北京大學(xué),2008.

[2] Dutton G.Improving locational specificity of map data——A multiresolution,metadata-driven app roach and notation[J].Int.J.Geographical Information Systems,1996,10(3）:253-268.

[3] GOODCH ILD M.Discrete global grids for digital earth[A].International Conference on Discrete Global Grids[C].California: Santa Barbara,2000.

[4] KIM ERL ING A,SAHR J.Comparing geometrical properties of global grids[J].Cartography and Geographic Information Science,1999,l26(4）:271-288.

[5] 鄧雪清.柵格型空間數(shù)據(jù)服務(wù)體系結(jié)構(gòu)與算法研究[D].中國人民解放軍信息工程大學(xué),2003.

[6] 杜瑩.全球多分辨率虛擬地形環(huán)境的金字塔模型研究[D].中國人民解放軍信息工程大學(xué),2006.

[7] 趙學(xué)勝.全球離散格網(wǎng)的空間數(shù)字建模[M].北京:測繪出版社,2007.

[8] 張永生.地球空間信息球面離散網(wǎng)格——理論、算法及應(yīng)用[M].北京:科學(xué)出版社,2007.

[9] 王家耀.地圖學(xué)原理與方法[M].北京:科學(xué)出版社,2007.

[10] 龔雪晶.基于Morton碼的圖像分裂合并算法研究[J].計算機工程與設(shè)計,2007,28(22）:5440-5443.

Global Subdivision Systems Based on Map Kilo-Grids

GENG Xiao-hui,CHENG Cheng-qi,SONG Shu-hua,L IDa-peng
(Institute of Remote Sensing and GIS,Peking University,Beijing 100871,China）

Kilo-grid is the grid fo rmed by the evenly spaced parallels of the p rojection coordinate axes.Re-p rojecting the kilogrids onto the earth ellip soid results in the subdivision of the earth ellipsoid,and this avoids the shape degeneration of the traditional geographic grid model in the high-latitudes and the irregularity p roblem of the meshes of the polyhedral-based discrete grid models effectively.Based on the quad-tree structure between the gridsof different levels,the acco rding mesh code is established on the basisof Mo rton Code and the transfo rmation between traditional geographic coordinates and the mesh codes and the calculation method of the adjacency relationship of meshes are imp lemented.In the end,related research is done on the distortion rules of Kilo-gridson the earth ellipsoidal surface.The results indicate that though there existsarea disto rtion w hen p rojecting themeshes to the ellipsoidal surface,the distortion is minor and follow s strict rules,w hich makes the Kilo-grid Model suitable as a basic Global Subdivision Model.

Global Subdivision System;Kilo-grids;Gause-Kruger p rojection;Morton code

P208

A

1672-0504(2010）02-0015-04

2009-10-27;

2009-12-22

國家重點基礎(chǔ)研究發(fā)展規(guī)劃(973）項目(61399）

耿曉暉(1979-）,男,碩士研究生,主要從事地理信息系統(tǒng)研究與工程應(yīng)用。E-mail:gengxiaohui2007@126.com